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Abstract 


Mathematical  modeling  of  biological  systems  is  crucial  to  effectively  and  efficiently  developing  treatments 
for  medical  conditions  that  plague  humanity.  Systems  of  differential  equations  are  the  traditional  tools  used 
to  theoretically  describe  the  spread  of  disease  within  the  body.  In  this  project  we  consider  the  dynamics  of 
the  Human  Immunodeficiency  Virus  (HIV)  in  vivo  during  the  initial  stages  of  infection.  Both  mathematical 
and  biological  results  support  the  idea  that  contact  with  the  HIV  retrovirus  does  not  automatically  imply 
permanent  infection.  Given  factors  such  as  the  CD4+  T-cell  growth  rate,  infection  rate,  and  viral  clearance 
rate,  it  is  possible  to  correctly  predict  the  end  viral  state  in  a  deterministic  model.  While  this  is  useful, 
such  a  model  lacks  the  randomness  inherent  in  physical  processes  and  parameter  estimation.  To  account  for 
this,  our  project  examines  both  discrete  and  continuous  stochastic  models  for  the  early  stages  of  HIV  infec¬ 
tion.  These  models  are  crafted  using  the  knowledge  of  biological  interactions  and  fundamental  mathematical 
principles.  We  also  examine  the  well-known  three-component  deterministic  model  in  greater  detail,  proving 
existence  and  uniqueness  of  the  solutions.  Furthermore,  we  prove  that  solutions  remain  biologically  meaning¬ 
ful,  i.e.,  are  positivity  preserving,  and  perform  a  thorough  stability  analysis  for  the  equilibrium  states  of  the 
system.  Finally,  we  develop  two  new  stochastic  models  and  obtain  extensive  numerical  results  to  measure  the 
probability  of  infection  given  the  transmission  of  the  virus  to  a  new  individual.  To  simulate  the  dynamics  of 
the  virus,  we  employ  a  number  of  computational  methods  for  ordinary  and  stochastic  differential  equations, 
including  Runge-Kutta  methods  and  the  Euler-Maruyama  scheme. 


Keywords:  mathematical  biology,  stochastic  differential  equations,  HIV,  differential  equations. 
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CHAPTER  1 

Introduction 


The  human  body  has  been  studied  for  hundreds  of  years,  and  the  knowledge  accumulated  to  date 
includes  the  intricate  molecular  level  mechanisms  that  determine  how  a  virus  infects  cells  and  reproduces. 
This  information  is  necessary  to  develop  treatments  for  maladies  ranging  from  Hepatitis  to  the  Human 
Immunodeficiency  Virus  (HIV)  and  Acquired  Immunodeficiency  Syndrome  (AIDS).  An  understanding  of  the 
accumulation  of  millions  of  these  interactions  across  a  timescale  of  medical  relevance  permits  the  creation 
of  better  treatment  methods.  As  recently  as  fifty  years  ago,  this  could  only  be  done  via  experiment  and 
repeated  trials.  Today,  through  the  advent  of  high  speed  computing,  viral  kinetics  can  be  simulated  using 
mathematical  models.  We  can  utilize  our  current  understanding  of  molecular  and  cellular  dynamics,  describe 
them  in  mathematical  terms,  and  simulate  interactions  within  the  body  many  thousands  of  times  to  examine 
the  course  of  a  virus  from  initial  infection  to  a  long-term  steady  state.  In  this  particular  project,  we  will 
derive  a  mathematical  model  of  HIV  infection  within  a  single  human  host,  examine  it,  simulate  it,  and 
improve  upon  it. 

The  general  principle  behind  any  sort  of  mathematical  model  is  to  examine  closely  the  interactions 
between  the  quantities  being  analyzed.  We  begin  with  a  small  number  of  initial  axioms,  and  then  follow  the 
implications  of  these  axioms  thoroughly  to  their  conclusion.  As  a  basis  for  our  model,  we  use  the  general 
biological  understanding  of  HIV  dynamics  that  is  widely-accepted  and  well  studied.  The  following  will  be  a 
brief  non-technical  discussion  of  the  populations  and  interactions  that  impact  the  modeling  of  HIV  within 
the  body. 

There  are  many  biological  operators  involved  in  the  interaction  between  HIV  and  cells  within  the  human 
body.  We  deal  with  a  small  portion  of  these  populations.  The  first  group  is  a  subset  of  the  population  of 
lymphocytes,  which  in  turn  are  a  type  of  white  blood  cell.  This  subset  is  known  as  CD4+  T-cells,  or  helper 
T-cells.  These  T-cells  have  a  variety  of  functions,  including  secreting  substances  that  stimulate  the  immune 
system,  acting  as  memory  agents  for  the  immune  system,  and  regulating  the  immune  response.  In  short, 
CD4+  T-cells  detect  and  direct  immune  system  responses  to  invading  bacteria  and  viruses.  Without  them, 
the  body  suffers  from  opportunistic  infections  that  are  greater  in  severity  and  duration  than  they  would  be 
if  the  CD4+  T-cells  were  otherwise  not  present.  HIV,  which  refers  in  this  case  to  the  virus,  not  the  disease 
or  symptoms  associated  with  it,  is  a  retrovirus  that  infects  helper  T-cells.  The  virus,  which  is  significantly 
smaller  than  the  T-cell  (120  nanometers  in  diameter  compared  with  7  micrometers  in  diameter),  breaches 
the  cell  wall  and  transports  its  RNA  into  the  T-cell  nucleus,  where  it  may  remain  dormant  for  a  time.  Upon 
activation,  the  T-cell  ceases  its  function  as  part  of  the  immune  system  and  instead  produces  additional  copies 
of  HIV.  These  infected  T-cells,  along  with  the  healthy  T-cells  and  free  floating  HIV,  are  the  populations 
with  which  we  are  concerned,  and  will  appear  within  our  mathematical  model. 

The  second  aspect  of  the  model’s  creation  involves  incorporating  the  interactions  between  these  popu¬ 
lations  and  deriving  their  corresponding  mathematical  representation.  For  each  population,  we  need  only 
concern  ourselves  with  the  ways  in  which  the  population  is  increased  and  how  members  of  this  group  are 
removed.  Healthy  T-cells  are  created  from  stem  cells  in  the  bone  marrow,  and  mature  in  the  thymus.  While 
production  of  T-cells  does  decrease  with  the  aging  of  the  human  body,  it  is  considered  constant  in  this 
project  for  two  reasons.  First,  there  is  no  method  other  than  this  production  that  can  possibly  affect  T-cell 
creation,  and  second,  the  time  scale  of  interest  within  the  model  is  sufficiently  small  to  consider  the  T-cell 
production  rate  as  constant.  When  considering  removal  of  these  cells,  we  note  that  T-cells  do  age  and,  in 
time,  die.  Within  the  model,  we  assume  that  each  T-cell  lives  for  roughly  the  same  amount  of  time,  and 


CHAPTER  1.  INTRODUCTION 


6 


thus,  the  death  rate  does  not  vary  over  the  entire  population.  Instead,  the  overall  number  of  T-cells  lost  in 
a  group  over  a  certain  period  of  time  is  proportional  to  the  number  of  T-cells  within  the  group.  The  other 
mechanism  through  which  the  population  of  healthy  T-cells  may  decrease  is  through  infection.  In  consid¬ 
ering  the  effects  of  infection,  we  utilize  an  interaction  term  that  arises  from  the  widely-used  “mass  action 
principle”  to  describe  the  transfer  which  occurs  between  populations  when  the  virus  infects  healthy  T-cells. 
This  mass  action  term  represents  the  idea  that  the  rate  of  interaction,  or  infection,  is  directly  proportional  to 
the  product  of  the  participating  populations,  namely  those  of  virions  (or  virus  particles)  and  healthy  T-cells. 
This  completes  the  interactions  for  the  healthy  T-cells.  We  note  that  the  only  way  to  add  to  the  population 
of  infected  T-cells  is  for  the  HIV  virus  particle  to  physically  infect  healthy  T-cells.  Beyond  this,  there  is 
no  other  way  of  creating  an  infected  T-cell.  Thus,  in  our  mathematical  model  we  can  use  the  same  mass 
action  term  to  describe  the  removal  of  the  healthy  T-cells  and  the  addition  of  infected  T-cells.  Similar  to 
the  healthy  T-cells,  infected  T-cells  die  off  or  are  cleared  by  the  immune  system  at  a  rate  proportional  to 
the  size  of  their  current  population.  The  virus,  while  produced  from  the  infected  T-cell,  does  not  cause  the 
destruction  of  the  infected  T-cell.  Thus,  we  do  not  need  to  include  such  a  transition  within  the  model. 

While  the  virus  production  rate  does  differ  from  cell  to  cell,  we  can  assume  that  the  aggregate  rate 
is  constantly  proportional  to  the  population  of  infected  T-cells.  This  is  the  only  mechanism  by  which  the 
virus  can  be  created.  Contrastingly,  there  are  two  manners  in  which  the  virus  can  be  removed,  or  cleared, 
from  the  body.  The  first  is  through  viral  infection  of  T-cells.  The  act  of  infecting  a  healthy  T-cell  must 
technically  remove  a  virus  particle  from  the  population  of  viruses  that  can  infect  further  T-cells.  However, 
when  considering  the  overall  population  quantities,  the  number  and  rate  of  viruses  lost  this  way  is  minute 
compared  to  other  methods  of  creation  and  destruction.  Hence,  we  will  omit  this  possibility  from  the  model. 
The  second  method  of  removal  is  known  as  viral  clearance.  It  is  the  removal  by  the  body  of  individual 
virus  particles,  and  is  performed  at  a  rate  proportional  to  the  current  amount  of  virus  particles  within  the 
body.  With  these  interactions,  we  can  construct  differential  equations  that  model  the  rates  of  change  of 
these  populations  by  utilizing  the  assumptions  discussed  above. 

For  all  of  the  interactions  we  have  described,  the  rates  of  change  were  proportional  to  a  population. 
In  the  equations  that  follow,  we  endow  each  of  these  interactions  with  a  proportionality  constant,  in  the 
necessary  units  to  guarantee  that  each  term  within  the  equation  is  measured  in  organisms  (cells  or  virions) 
per  unit  time.  For  instance,  the  creation  rate  of  the  healthy  T-cells  will  be  represented  by  the  constant 
A  >  0,  whose  unit  of  measurement  is  T-cells  per  unit  time.  Analogously,  the  rate  coefficient  for  the  creation 
of  infected  T-cells  is  measured  in  unites  that  are  the  inverse  of  virus  particles  multiplied  by  a  unit  of  time. 
Thus,  when  we  sum  up  the  interactions  for  each  population,  we  arrive  at  a  growth  rate.  Mathematically, 
the  rate  of  change  of  a  quantity  is  given  by  its  derivative,  and  hence  a  differential  equation  results.  In  our 
case,  the  system  is  known  as  the  Three  Component  Model  (3CM)  16  . 

Before  stating  the  equations,  we  attach  variable  names  to  each  of  the  populations  under  consideration. 
In  the  notation  used,  the  population  of  healthy  T-cells  is  denoted  by  T,  the  infected  T-cells  by  /,  the  virus 
population  by  V,  and  the  derivative  of  any  population  X  with  respect  to  time  is  indicated  by  .  as  is 
standard.  There  are  six  different  proportionality  constants,  given  in  the  table  below. 


Coefficient 

Relevant  Biological  Process 

Mean  Value 

Units 

A 

Healthy  T-cell  growth 

0.1089 

cells 

day 

Healthy  T-cell  death 

0.01089 

1 

day 

k 

Infection 

1.179  xlO-3 

1 

virions  X  day 

S 

Infected  T-cell  death 

.3660 

1 

day 

P 

Virus  production 

1426.8 

virons 

cellsxday 

c 

Viral  clearance 

3.074 

1 

day 

We  note  several  important  characteristics  of  these  coefficients.  Since  the  coefficients  are  taken  from 
biological  considerations,  we  assume  that  they  are  all  nonnegative,  and  we  fix  their  units  such  that  the 
end  result  of  each  equation  is  ultimately  represented  in  units  of  or  V*T ™es .  The  units  for  each  T-cell 
population  is  simply  “cells”,  while  the  units  for  the  virus  population  is  “virons”. 
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Figure  1.1.  Illustration  of  3CM 

When  we  combine  the  multiple  interactions  that  govern  the  change  in  population  of  the  healthy  T-cells, 
including  their  natural  growth  rate,  natural  death  rate,  and  infection  rate,  we  arrive  at  the  differential 
equation 

(1)  f^X-.T-kTV. 

This  equation  mathematically  describes  the  biological  interactions  presented  above.  A  similar  equation  is 
derived  for  the  infected  T-cells,  namely 

(2)  ft=kTV~SI' 

Here,  it  is  important  to  note  that  the  “ kTV ”  terms  in  the  first  and  second  equations  represent  the  same 
interaction  and  yield  the  same  number  of  cells  transferred  from  the  healthy  T-cell  population  to  the  infected 
T-cell  population.  This  occurs  because  the  transfer  of  cells  between  these  populations  occurs  at  a  one  to  one 
ratio. 

Next,  we  consider  the  rate  of  change  of  the  virus  population  and  derive  the  corresponding  differential 
equation  to  represent  its  dynamics: 

(3)  ^=U-cK 

The  diagram  ([l])  gives  a  visual  representation  of  the  differential  equations.  We  see  each  of  the  populations 
within  a  circle,  while  the  arrows  running  to  and  from  each  circle  describe  how  they  interact  and  change. 


Of  course,  the  unknown  populations  within  these  equations  depend  upon  one  another,  and  hence  they 
are  not  examined  independently.  Just  as  the  interaction  of  populations  within  the  body  serves  to  connect  the 
behavior  of  virions  and  T-cells,  the  functions  within  the  differential  equations  above  are  also  coupled.  They 
describe,  for  instance,  the  notion  that  an  increase  within  the  virus  population  will  cause  a  similar  increase 
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within  the  infected  T-cell  population.  Because  of  these  interactions,  we  study  the  equations  together  and 
look  at  properties  of  the  system  of  three  equations,  rather  than  an  isolated  equation.  Thus,  the  equations 
together  are  known  as  the  “Three  Component  Model.” 


(3CM) 


dT 

dt 

dl 

— 

dt 
dV 
'  dt 


A  -fj,T-  kTV 
kTV  -  SI 
pi  —  cV. 
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CHAPTER  2 

The  Deterministic  Three  Component  Model 

2.1.  Existence  and  Uniqueness  of  Solutions 

The  most  fundamental  step  in  examining  the  Three  Component  Model  is  to  show  that  a  solution  to  the 
initial-value  problem  does,  in  fact,  exist,  and  that  this  solution  is  unique.  To  prove  this,  we  utilize  the  classic 
theorem  of  Picard  and  Lindelof,  which  is  stated  in  3 1 . 

Theorem  2.1.1  (Picard-Lindelof).  Let  n  £  N  and  yo  £  R"  be  given.  Assume  the  function  f  :  Rn  xl-} 
R"  is  locally  Lipschitz  in  its  first  argument  and  continuous  in  its  second  argument.  Then  there  exist  t*  >  0 
and  a  unique  function  y  :  [0,  f*]  — >  R71  satisfying 

y'(t )  =  f{y(t),t ) 

for  every  t  £  [0,  t*]  and  the  initial  condition  y( 0)  =  yo. 

In  the  case  of  the  three  component  model,  we  have: 


(4) 

y  = 

T 

I 

and  the  autonomous  system  f(y)  = 

A  —  pT  —  kTV 
kTV  -  SI 

V 

pi  —  cV 

We  note  that  since  /  has  a  continuous,  bounded  derivative  on  any  compact  subset  of  R3  that  /  is  locally 
Lipschitz  in  y.  Hence,  by  the  Picard-  Lindelof  Theorem,  there  exists  a  unique  solution,  y(t),  to  the  ordinary 
differential  equation  y'(t)  =  f(y(t),t)  on  some  time  interval  [0,  f*]. 


2.2.  Gronwall’s  Inequality  and  Applications  to  the  Three  Component  Model 

Gronwall’s  inequality  is  a  tool  of  mathematical  analysis  that  allows  us  to  bound  solutions  of  differential 
equations.  In  many  cases,  we  cannot  explicitly  obtain  formulas  for  solutions,  so  it  must  suffice  to  derive 
bounds  on  their  growth.  Gronwall’s  Inequality  permits  us  to  do  just  this. 

2.2.1.  Gronwall’s  Inequality. 

Theorem  2.2.1.  Let  [a,  b]  be  an  interval  and  f  and  g  be  continuous  functions  on  [a,  b]  with  f  differentiable 
on  [a,  b}.  If  f  satisfies  the  differential  inequality: 

(5)  ^  =  f'(t )  <  f{t)g(t) 

for  every  t  £  [ a,b ]  then 

(6)  f(t)<f(a)eK^ds 

for  every  t  £  [ a,b ]. 

The  proof  to  this  theorem  is  well  known,  and  can  be  found  in  6  . 
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2.2.2.  Application  to  the  Three  Component  Model. 

Theorem  2.2.2  (Boundedness  and  Positivity  for  the  Deterministic  Three  Component  Model).  Let  to  >  0 
be  given.  In  the  three  component  model,  if  the  initial  conditions  satisfy  T( 0)  >  0,  1(0)  >  0,  and  U(0)  >  0, 
then  for  all  t  €  [0,  fo],  T(f), /(f),  and  V(t)  will  be  bounded  and  remain  positive. 

Proof.  We  assume  that  T(t),  I(t),  and  V(t)  initially  have  positive  values.  From  the  previous  theorem, 
there  exists  a  t*  such  that  the  solution  exists  on  [0,f*].Let  us  denote  by  T*  the  largest  time  t  for  which 
T(t)  >  0,  /(f)  >  0,  and  V(t)  >  0,  or  more  precisely 

T*  =  sup{t  G  [0,t*]  :  T(s),I(s),V(s)  >  0  for  all  s  G  [0,t]}. 

Then  on  the  interval  [0,T*]  we  can  make  estimations  stemming  from  Gronwall’s  inequality. 

Recall  that  all  constants  in  the  equations  are  positive.  Using  Gronwall’s  inequality,  we  can  place  lower 
bounds  on  jjj  and 


and 


^  =  kTV  -  61  >  -SI 
dt 

dV 

— —  =  pi  —  cV  >  —cV 
dt 


Using  an  integrating  factor,  we  rewrite  these  differential  inequalities  to  find 

/(f)  >  e~5t  >  0 

for  f  €  [0,  T*]  and 


V(t)> 


>  0 


for  f  G  [0,  T*].  Similarly,  we  can  place  an  upper  bound  on  ^ : 


Solving  for  T  yields: 


lrri 

—  =  A  -  pT  -  kTV  <  A. 
dt 


T(f)  <T(0)  +  Af  <ci(l+f). 


Define  the  constant  ci  such  that  ci  >  max(A,T(0)).  We  can  sum  the  equations  for  dd.  and  'jV  and  place 


bounds  on  this  sum  so  that 


—  (/  +  V)  =  kTV  -  61  +  cV  -  pi  <  kTV  +  pi. 
dt 


Recall  that  we  have  a  bound  on  T,  so  we  can  substitute: 

d 


dt 


(I+V)  <  kci(l  +  t)V  +pl  <  c2(l  +  f)(/  +  V) 


where  C2  >  max(/cci ,  p) .  Solving  the  differential  equation  yields 

(/  +  V)(t)  <  (1(0)  +  V(0)  +  c2)e*2  <  c3e‘2 

for  f  G  [0,  T*]  where  c3  >  3  max((/(0),  V(0),  c2).  Since  I(t)  is  positive,  we  can  place  an  upper  bound  on  V: 

c3e‘2  >  (/  +  U)(f)  >  V(t) 

We  also  know  that  since  V(t)  is  positive,  so  is  /(f). 

c3e<2  >  (I  +  V)(t)  >  /(f) 

With  these  bounds  in  place,  we  can  now  examine  T(t): 

dT  2 

—  =  A  -pT-  kTV  >  -pT  -  kTV  >  -pT  -  kc3ef  T 

dt 

>  -c4(l  +  e‘2)T 
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for  t  G  [0,  T*],  where  C4  >  max(/z,  kc3).  Shifting  that  last  term  to  the  other  side  of  the  equation  yields 

^+c4(l  +  e‘2)T>0 

Since  we  know 

^(T(t)+eC4/o(1+er2^))>0< 

at 

then  we  find 

(7)  T(f)  >  T(0)e_C4  /o  (1+eT2dr)  >  0 
for  t  G  [0,  T*]. 

Thus,  the  values  of  T,  I,  and  V  stay  strictly  positive  for  all  of  [0,T*],  including  the  time  T*.  By 
continuity,  there  must  exist  a  t  >  T*  such  that  T(t),  I(t),  and  V(t)  are  still  positive.  This  leads  to  a 
contradiction,  and  shows  that  T,  J,  and  V  are  strictly  positive  on  the  entire  interval  [0,T*].  Additionally,  on 
this  same  interval,  all  of  the  functions  remain  bounded,  so  the  interval  of  existence  can  be  extended  further. 
In  fact,  the  bounds  on  T,  /,  and  V  derived  above  hold  on  any  compact  time  interval.  Thus,  we  may  extend 
the  time  interval  on  which  the  solution  exists  to  [0,  to]  for  any  to  >  0  and  from  the  above  argument,  the 
solutions  remain  both  bounded  and  positive  on  [0, to].  □ 

With  this  knowledge,  we  remain  secure  in  the  soundness  of  our  model  and  know  with  certainty  from  a 
mathematical  perspective  that  it  retains  biological  validity.  This  also  shows  that  once  infected,  it  is  entirely 
possible  that  the  virus  may  continue  to  exist  beneath  a  detectable  threshold  without  doing  any  damage. 

2.3.  Stability  and  Equilibrium  of  the  Deterministic  Three  Component  Model 

2.3.1.  Equilibria  of  the  System.  In  order  to  fully  understand  the  three  component  model,  it  is 
necessary  to  first  study  the  equilibrium  points. 

Definition  2.3.1.  Consider  the  differential  equation  y'(t)  =  f(y(t),t).  A  point  is  an  equilibrium  point 
=  f{y{t),t )  =  0  for  all  t  G  R. 

In  our  case,  an  equilibrium  point  is  a  set  of  three  values  corresponding  to  (T,  /,  V)  with  the  property 
that  if  the  system  begins  at  these  values,  it  will  remain  there  indefinitely.  In  other  words,  the  populations 
are  unchanging  from  those  values,  so  the  rate  of  change  for  each  population  is  zero.  Knowing  that  is  zero 
is  what  allows  us  to  solve  for  the  precise  values.  In  the  case  of  the  Three  Component  Model,  there  exist 
exactly  two  equilibrium  points.  We  discover  them  by  setting  dJf ,  dJ ,  and  ^  equal  to  zero  and  solving  the 
resulting  equations  for  T,  /,  and  V.  From  a  biological  perspective,  we  can  categorize  these  points  to  be  when 
the  HIV  virus  is  either  extinct  from  the  body,  i.e. ,  I  =  V  =  0,  or  when  the  virus  persists  within  the  body 
(I  0,  V  ^  0)  as  t  grows  large. 

2.3.2.  Determining  Equilibrium  Values.  The  first  such  equilibrium  exists  when  I  =  0  and  V  =  0. 
This  results  in  a  simplification  of  the  above  system  of  equations  to: 

(8)  0  =  A  -  fj,T, 

which  in  turn  shows  that  the  ordered  triplet  (T,  /,  V)  =  (^,0,0)  is  an  equilibrium  solution.  This  particular 
equilibrium  point  is  also  known  as  viral  extinction,  since  there  are  no  virus  particles  or  infected  cells.  We 
will  refer  to  this  point  as  P\. 

To  find  the  second  equilibrium,  we  assume  I  ^  0  and  V  0.  After  setting  AL  =  -L  =  dv.  =  q,  solving  for 
/  in  the  second  equation  reveals  that  I  =  while  solving  for  I  in  the  third  equation  gives  I  =  ^ ,  It  follows 
that  since  V  ^  0,  T  =  ||.  Substituting  this  value  of  T  into  the  first  equation  yields  V  =  and  further 

substitution  shows  I  =  |  Thus,  there  is  a  second  equilibrium  at  (T,I,V)  =  (||,  |  ^  ^ ).  Since 

there  are  distinct  presences  of  virus  particles  and  infected  T-cells,  we  refer  to  this  point  as  viral  persistence 
and  abbreviate  the  point  as  P2  ■  In  terms  of  biology,  we  can  say  P\  is  the  case  in  which  an  infection  exists  for 
a  short  period  of  time,  then  is  removed  from  the  body  by  natural  means.  The  virus  does  not  persist.  The 


2.3  Stability  and  Equilibrium  of  the  Deterministic  Three  Component  Model 


12 


second  case,  where  the  system  of  equations  tends  to  P2,  denotes  that  situation  where  the  body  is  unable  to 
clear  the  infection  by  itself.  If  this  ends  up  being  the  case,  than  after  a  certain  period  of  time,  the  Three 
Component  Model  loses  its  applicability  as  the  infection  takes  a  deeper  hold  on  the  body.  More  complex 
models,  which  consider  latent  infection,  effects  of  macrophages,  or  cytotoxic  immune  response,  are  then 
required  to  describe  the  spread  of  HIV  within  the  body  and  its  development  towards  AIDS. 

If  the  system  takes  on  the  values  of  Pi  or  P2  at  any  time,  it  will  remain  at  those  points  indefinitely. 
Unless  the  initial  conditions  for  the  system  are  set  at  an  equilibrium,  in  general,  the  system  may  not  actually 
obtain  these  values.  Instead,  the  system  will  likely  approach  them  asymptotically,  tend  to  infinity,  or  cycle 
between  certain  values.  In  the  case  of  the  Three  Component  Model,  we  will  see  that  the  system  approaches 
each  equilibrium  asymptotically.  This  is  a  concept  known  as  asymptotic  stability: 

Definition  2.3.2.  An  equilibrium  solution  y*  to  an  ordinary  differential  equation 

y'(t )  = 

is  said  to  be  (locally)  asymptotically  stable  if 

lim  |  y(t)  -  y*\  =  0 

t—>  OO 

whenever  the  initial  values  yo  are  sufficiently  close  to  y* . 


For  linear  ODEs,  it  is  well-known  that  the  stability  properties  depend  only  upon  the  eigenvalues  of  the 
system.  However,  our  model  (3CM)  is  nonlinear,  and  thus  we  must  rely  on  linearization  and  a  theorem  of 
Hartman  &  Grobman  7  to  unify  the  behavior  of  the  linear  and  nonlinear  systems. 


2.3.3.  Linearization  and  the  Jacobian  of  the  System.  Linearization  is  a  key  concept  when  exam¬ 
ining  the  equilibrium  stability  of  a  system  of  differential  equations.  In  this  case,  we  approximate  the  system 
of  differential  equations  with  a  linear  system  at  the  points  P\  and  P2.  This  stability  analysis  perturbs  the 
system  from  equilibrium  and  examines  if  the  system  returns  to  the  original  equilibrium.  In  order  to  linearize 
the  system  we  need  to  first  compute  the  Jacobian  matrix  for  the  system.  If  a  system,  5,  consists  of  n 
functions  of  m  variables,  i.e.,  F\(x\, ...,  xm), ...,  Fn(x i, ... xm )  then  the  Jacobian  matrix,  Jg,  is  a  matrix  of  the 
partial  derivative  of  each  function  with  respect  to  each  variable: 


rdF\ 

d-Fi  -| 

dxi 

‘  dXm 

dFn 

dFn 

-  dxi 

dxm  - 

The  matrix  J$  provides  a  linear  approximation  of  a  system  at  any  given  point,  and  when  evaluated  at 
an  equilibrium  point  P,  Js(P)  also  encodes  information  above  the  nonlinear  system. 

It  is  necessary  to  evaluate  a  Jacobian  matrix  at  both  Pi  and  P2  and  examine  its  corresponding  eigen¬ 
values,  since  analysis  of  the  eigenvalues  of  the  Jacobian  matrix  evaluated  at  a  equilibrium  gives  insight  into 
the  stability  properties  of  that  equilibrium.  Specifically,  if  the  eigenvalues  of  the  Jacobian  matrix  at  a  point 
have  negative  real  parts,  then  that  point  is  classified  as  an  attractor,  meaning  small  perturbations  from  the 
equilibrium  result  in  the  system  returning  to  that  point  over  time.  On  the  other  hand,  if  one  or  more  of 
the  eigenvalues  have  positive  real  part,  then  small  perturbations  from  equilibrium  result  in  magnifications 
of  those  disturbances,  and  the  system  shifting  away  from  the  point.  This  point  would  then  be  known  as  a 
“repeller.”  For  a  more  detailed  look  at  eigenvalues  and  their  effects  on  systems  of  differential  equations,  see 


In  the  case  of  (3CM),  the  Jacobian  matrix  is  given  by 


J-scm 


-kV  —  y  0  -fcT 
kV  - 6  kT 
0  p  —  c 


It  is  possible  to  determine  the  signs  of  the  eigenvalues  of  the  Jacobian  using  a  theorem  of  Routh  and  Hurwitz 
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2.3.4.  Routh-Hurwitz  Theorem. 

Theorem  2.3.3  (Routh-Hurwitz  Criteria  for  a  Characteristic  Polynomial).  Given  the  polynomial  P(x)  = 
xn  +  a\x n~1  +  •  •  •  +  an-\X  +  an  where  each  ai  is  a  constant  real  coefficient,  define  the  n  Hurwitz  matrices 
using  the  coefficients  a*: 


„  1 

(Z3  1 

o' 

Hi  —  [03]  ,  H2 

CL\  1 

O3  0,2 

,h3  = 

a3  a2 

a  3 

CZ5  04 

03. 

ai  1 

0  0 

. .  O' 

03  «2  ai  1 

. .  0 

05  CI4  CI3  02 

. .  0 

Hn  = 

<27  ae  a<. 5  04 

. .  0 

_0 

0 

0  0 

•  •  &n_ 

All  of  the  roots  of  P(x)  have  negative  real  part  if  and  only  if  the  determinants  of  all  of  the  Hurwitz  matrices 
are  positive. 

If  we  consider  the  case  n  =  3,  this  theorem  simplifies  significantly:  Given  the  polynomial  x 3  +  a\x2  + 
a2X  +  03,  the  three  Hurwitz  matrices  are  as  follows: 


_ 

ai 

1 

O' 

a  1 

1 

:  H:i  = 

a3 

0,2 

Oi 

03 

a2_ 

0 

0 

o3_ 

If  the  determinants  of  Hi,  H2 ,  and  H3  are  all  positive,  then  the  roots  of  the  characteristic  equation,  i.e., 
the  eigenvalues,  are  all  negative  or  have  negative  real  parts.  We  can,  in  fact,  make  further  simplifications. 


Corollary  2.3.4.  Based  on  the  Routh-Hurwitz  criteria,  if  a\  >  0,aiO2  >  03,  and  03  >  0  then  all  of 
the  eigenvalues  have  negative  real  part,  implying  that  this  particular  equilibrium  is  asymptotically  stable. 

PROOF.  Since  Det(Hi)  =  a±,  Det(H2)  =  0302  —  03,  and  Det(fl3)  =  030203  —  afi1 ,  it  is  equivalent  to  the 
Routh-Hurwitz  criteria  to  say  that: 

(1)  oi  >  0,  so  that  Det(Ri)  is  positive 

(2)  0102  >  03,  so  that  Det(I?2)  is  positive,  and 

(3)  a3  >  0,  so  that  Det(R3)is  positive. 

This  last  statement  follows  from  the  first  two.  Since  0302  >  03,  and  Det(H3)  =  030203  —  a2  =  03(0302  — 
03),  we  see  that  03(0302  —  03)  is  positive  if  and  only  if  03  >  0. 

□ 


As  we  will  see  below,  the  signs  of  the  roots  of  this  polynomial  depend  only  on  a  single  number  known 
as  the  reproductive  constant,  R.  Therefore,  it  is  sufficient  to  examine  R  =  and  the  effect  it  has  on 
the  coefficients  of  the  characteristic  equation  to  determine  how  R  influences  stability  of  equilibrium.  This 
means  that  in  a  deterministic  model,  we  can  simply  examine  the  value  of  this  coefficient  to  determine  viral 
persistence  or  extinction.  This  is  a  remarkable  result,  and  estimates  of  the  persistence  of  HIV  upon  initial 
infection  have  been  made  solely  by  Monte-Carlo  simulations  generating  different  R-values 
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2.3.5.  The  Relation  between  R  and  Equilibrium  Stability.  Here,  we  state  two  theorems  that 
demonstrate  the  relationship  between  the  R-value  and  which  equilibrium  is  asymptotically  stable. 

Theorem  2.3.5.  P3  is  an  asymptotically  stable  equilibrium  if  and  only  if  R.  is  less  than  1. 
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PROOF.  Evaluating  the  Jacobian  at  Pi  =  (^,0,0)  results  in 


J3Gm{Pi)  = 


-a*  0  -T 


-8  ^ 

p  —C  _ 


The  characteristic  equation  of  Jzcm(P^)  is: 


(— p  —  x )  (—8  —  x)(—c  —  x) 


P  / 


=  0 


which  expands  to: 
where 


x  T  d\X  T  d^x  T  03  —  0 

a±  =  c  +  5  +  fxa.2  =  c8  +  cp  +  8p  — -^—03  =  cSp  —  kpX. 


If  the  eigenvalues  have  negative  real  part,  then  the  Routh-Hurwitz  criteria  are  satisfied,  which  implies 
that  03  >  0.  This  means  c8p  —  kpX  >  0  and  it  follows  that  c8p  >  kpX  leading  to  1  >  =  R.  Simply  put, 

an  asymptotically  stable  equilibrium  at  Pi  implies  that  R  <  1. 

Conversely,  if  R  <  1,  then  03  >  0,  and  since  all  the  coefficients  are  positive,  it  is  clear  that  oi  >  0. 
Additionally,  since  cSp  >  kpX,  then  c8  >  and  a 2  =  cS  +  cp  +  Sp - ^  >  cS  +  cp  +  Sp  —  cS  =  cp  +  cS.  So, 


0,102  =  (c  +  5  +  p)  [  cS  +  cp  +  8  p  — 


kpX 


>  (c  +  8  +  p){cp  +  8p)  >  cSp  >  cSp  —  kpX  =  03 


that  is,  R  <  1  implies  Pi  is  an  asymptotically  stable  equilibrium  for  (3CM). 


□ 


Thus,  “R  <  1”  and  “Pi  is  an  asymptotically  stable  attracting  equilibrium”  are  equivalent  statements 

Theorem  2.3.6.  P2  is  an  asymptotically  stable  equilibrium  if  and  only  if  R  is  greater  than  one. 

Proof.  The  analysis  for  P2  =  (||,  y  —  fy,  yy  —  y)  is  similar  to  that  of  P\.  The  analysis  begins  by 
linearizing  the  system  (3CM|  about  P2  and  examining  the  characteristic  equation.  The  Jacobian  is  slightly 
more  complicated  in  this  case: 


0 

-Hgr 

JsCm{P2)  = 

-5 

0 

p 

— C 

■-O  0 

-(f)' 

kXP  n  A 

cd 

c5  k  ° 

P 

0  p 

—c 

and  results  in  the  characteristic  equation 
’  — kXp 


cS 


-  x)(—c  -x)-  c8)^  +  -P-  -  pSjp  =  0 


with  expanded  form 
where 


x3  +  a\x2  +  a-2.x  +  03  =  0 


r  kXp  kXp  kXp  ,  r 

oi  =  c  +  8  H - —  02  =  — — I - 03  =  kXp  —  cdp. 

cd  8  c 


Since  the  Routh-Hurwitz  criteria  apply  in  this  case  as  well,  it  is  sufficient  to  show  that  R  >  1  implies 
oi  >  0,  0102  >  a 3,  and  03  >  0  and  vice  versa. 
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It  is  clear  that  if  the  Routh-Hurwitz  criteria  are  satisfied  then  R  >  1,  since  CI3  =  kXp  —  c5p  >  0  43-  kXp  > 
cdp  <t=>  =  R  >  1.  Additionally,  since  all  the  coefficients  in  the  system  are  positive,  cti  will  be  positive  as 

well.  Finally, 


aia2  = 


(c  +  <5  + 


>  kpX  >  kpX  —  c5p  =  a 3. 


□ 


Thus,  R  >  1  implies  that  P2  is  an  asymptotically  stable  equilibrium,  and  in  fact,  is  both  a  necessary 
and  sufficient  condition.  This  analysis  reveals  one  very  important  fact  about  the  overall  system:  the  long 
term  behavior  of  the  system  depends  only  on  the  value  of  R.  If  R  is  greater  than  one,  then  the  system  tends 
towards  an  end  state  with  a  non-zero  population  of  infected  cells  and  virions  (persistence),  but  if  R  is  less 
than  one,  then  the  final  equilibrium  is  a  state  with  no  virus  or  infection  (extinction).  Finally,  we  remark 
that  in  the  rare  event  that  R  =  1,  the  equilibria  Pi  and  P2  are  identical.  From  the  analysis  above,  then,  it 
follows  that  this  joint  equilibrium  is  also  locally  asymptotically  stable. 

The  figure  (2.1)  serves  as  an  example  that  illustrates  what  the  equilibriums  looks  like  when  R  changes. 
The  top  graph  shows  what  happens  when  R  is  less  than  one,  while  the  lower  graph  shows  what  happens  when 
R  is  greater  than  one.  While  the  coefficients  are  not  biologically  relevant,  the  figure  remains  mathematically 
accurate. 


Viral  Extinction  A 


Figure  2.1.  Graphs  of  the  3CM  with  illustrative  coefficients  for  R  <  l(top)  and  R  >  1  (bottom). 
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2.4.  Numerical  Results 

The  simulation  of  the  Three  Component  Model  was  developed  in  MATLAB,  and  there  were  several 
different  paths  taken  to  get  final  probability  distributions,  and  each  of  them  differ  in  their  method  of  com¬ 
putation. 

2.4.1.  General  Method  of  Coding.  All  of  the  code  used  to  simulate  the  Three  Component  Model 
shares  the  same  basic  characteristics.  In  each  individual  trial,  initial  conditions  and  coefficients  are  chosen, 
the  R-value  is  calculated,  then  the  system  of  equations  is  simulated.  This  is  done  deterministically,  discretely, 
or  stochastically.  The  end  states  of  each  simulation  are  recorded  and  compared  with  the  viral  detection 
threshold,  and  in  some  cases,  we  determine  to  what  degree,  if  any,  correlation  exists  between  the  end  state, 
the  R-value,  and  the  initial  conditions.  We  also  examine  what  the  average  sample  path  and  average  end 
states  look  like. 


Figure  2.2.  A  flow  chart  representation  of  the  functions  involved  in  coding 


2.4.2.  Common  Features  of  Simulations.  There  are  several  functions  used  in  MATLAB  that  are 
common  to  nearly  every  simulation.  They  all  stem  from  a  specially  designed  random  number  generator. 
This  function  permits  the  generation  of  normal  distributions  that  are  truncated  on  one  or  both  ends.  Known 
as  truncatednormal  .m  throughout  the  project,  we  use  this  function  to  create  positive  or  non-negative 
coefficients  and  initial  conditions  that  are  both  statistically  representative  and  biologically  meaningful.  There 
are  child  functions  of  this  random  number  generator  that  create  truncated  normally  distributed  starting 
conditions  and  coefficients. 

In  addition  to  functions  used  throughout  all  the  simulations,  several  numerical  values  remain  the  same 
regardless  of  the  method  of  simulation  used.  For  instance,  there  exists  a  lowest  possible  amount  of  virus 
detectable  by  medical  technology.  40-75  virus  particles  per  micro-liter  of  blood  is  the  range  of  values  given 
by  scientific  literature  We  use  this  lower  value,  40  Partl^les  as  0ur  threshold.  We  also  run  each  simulation 
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from  initial  infection  to  sixty  days  beyond  that  point.  At  that  point,  the  populations  are  stable,  and  a 
determination  of  viral  extinction  can  be  made. 

Coding  the  three  component  model  deterministically,  while  not  fully  representative  of  the  randomness 
inherent  in  real  life,  still  offers  two  distinct  advantages.  First,  this  type  of  simulation  is  straightforward  from 
a  computational  perspective,  and  second,  this  simulation  allows  us  to  isolate  the  effects  of  coefficients  and 
initial  conditions  to  a  much  greater  degree  than  permissible  in  the  stochastic  simulations.  Since  we  recall  that 
the  equations  for  the  deterministic  model  and  for  the  expected  value  of  the  stochastic  model  are  the  same, 
we  can  imagine  this  simulation  as  looking  at  the  average  of  an  infinite  number  of  stochastic  simulations. 

2.4.3.  The  Basic  Deterministic  Model.  The  Deterministic  Model,  when  simulated  in  its  most  basic 
sense,  is  simulated  over  one  set  of  coefficients  and  initial  conditions.  The  system  of  differential  equations, 
along  with  the  initial  values  and  time  interval  over  which  it  is  to  be  simulated  are  given  to  the  MATLAB 
function  ode45 .  m,  which  uses  a  fifth  order  Runge-Kutta  method  to  generate  the  sample  path  for  each 
population.  With  one  trial,  it  is  a  simple  matter  to  compare  the  end  virus  population  to  the  detection 
threshold  and  determine  persistence  or  extinction  from  that. 


■e 

« 

£  o 


200  400 

Time  (days) 


200  400 

Time  (days) 


Figure  2.3.  Two  simulations  of  the  deterministic  3CM  with  mean  standard  ICs 


Figure [23] shows  what  one  run  of  the  deterministic  code  looks  like.  The  red  line  represents  the  detectable 
virus  threshold.  Note  that  this  demonstrates  the  two  possible  end  states  for  the  system.  We  can  either  have 
viral  persistence,  as  demonstrated  by  the  lower  set  of  figures,  or  we  can  have  viral  extinction,  as  demonstrated 
by  the  upper  set  of  figures.  These  figures  used  the  values  of  the  coefficients  that  would  maximize  and  minimize 
R,  in  order  that  a  clear  distinction  could  be  made  about  how  these  coefficients  affect  the  end  results. 

2.4.4.  Deterministic  Model  with  Random  Valued  Coefficients.  The  next  step  in  modeling  the 
deterministic  system  is  to  randomize  the  coefficients  used  in  the  equations.  We  recall  that  these  coefficients 
vary  from  person  to  person,  and  are  representative  of  various  biological  traits.  Because  of  this,  certain 
restrictions  and  conditions  are  placed  upon  them.  First,  we  know  that  in  order  to  be  biologically  meaningful, 
we  must  have  these  values  be  positive.  It  would  make  no  sense,  for  instance,  for  the  T-cell  growth  rate  to  be 
inherently  negative.  Second,  we  assume  these  values  to  be  normally  distributed.  Therefore,  when  we  draw 
upon  these  coefficients,  we  take  the  mean  and  standard  deviation  given  by  existing  data,  and  use  a  function 
in  MATLAB  designed  to  draw  a  random  value  from  a  truncated  normal  distribution. 
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For  each  of  10000  trials,  we  generate  a  set  of  coefficients  from  these  truncated  normal  distributions.  With 
these  values,  the  standard  initial  conditions,  and  the  standard  time  interval,  we  use  the  built  in  MATLAB 
function  ode45 .  m  again  to  computationally  solve  the  deterministic  system.  We  record  the  sample  paths  for 
each  run,  and  those  with  R- values  greater  than  one  are  colored  in  cyan  (or  light  blue),  while  those  sample 
paths  whose  equations  have  R- values  less  than  one  are  colored  magenta  (or  purple). 


Figure  2.4.  1000  simulations  of  the  deterministic  model  with  different  coefficients.  Stan¬ 
dard  IC 


In  the  end,  the  code  gives  an  array  of  data  consisting  ofthe  end  population  values  and  the  R-value  for  each 
trial  .  We  place  the  data  into  histograms  to  give  a  better  idea  of  how  many  end  states  reach  which  value. 
Again,  when  we  look  at  the  viral  end  states,  the  red  line  represents  the  detection  threshold.  In  this  particular 
example,  we  see  that  we  have  over  2000  cases  (2140  to  be  exact)  in  which  the  virus  population  is  undetectable 
after  sixty  days. 


End  Population  Value  End  Population  Value  End  Population  Value 

Figure  2.5.  1000  simulations  of  the  deterministic  model  with  different  coefficients.  Stan¬ 
dard  IC 


Of  the  coefficient  sets  generated,  211  had  corresponding  R- values  less  than  one.  This  is  similar  to 

whose  estimation  for  R-values  of  less  than  one  ranged  around  2%. 


results  of  Tuckwell  and  Shipman  17 


all  of  these  sets,  the  end  virus  population  was  beneath  the  detectable  threshold. 
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CHAPTER  3 


Stochastic  Differential  Equations  and  Ito’s  Lemma 


3.1.  The  Birth  Death  Model 


While  knowledge  of  the  deterministic  Three  Component  Model  is  useful,  our  goal  is  to  extend  this  to  a 
stochastic  model.  In  order  to  construct  the  stochastic  three  component  model,  we  first  create  a  stochastic 
model  of  a  simpler  system.  This  will  serve  to  illustrate  some  ideas  and  demonstrate  that  they  are  special 
cases  of  general  theorems  we  use  later  to  create  a  stochastic  model  of  the  Three  Component  Model. 

A  stochastic  system  better  represents  the  randomness  inherent  in  everyday  life.  When  used  with  one 
set  of  coefficients,  the  deterministic  model  tells  us  what  happens  to  one  patient  one  time  after  an  initial 
infection.  On  the  other  hand,  with  a  stochastic  model,  we  see  the  probability  distribution  of  this  one  person’s 
viral  state.  This  is  a  much  more  valuable  piece  of  information,  because  from  it,  we  can  see  not  only  which 
end  states  are  more  likely,  but  how  likely  they  are. 

Let’s  begin  by  considering  the  ordinary  differential  equation: 


(9) 


dX 

dt 


=  bX  -  dX 


X(0)  =  x0 


where  b,  d  >  0  and  Xq  >  0.  This  is  known  as  the  birth-death  model. 


3.1.1.  A  random  walk  related  to  the  birth-death  model.  We  wish  to  create  a  random  walk  that 
behaves  in  such  a  way  that  the  expected  value  of  the  random  walk  at  time  t  is  related  to  the  differential 
equation  (|9|.  To  do  this,  we  first  define  a  discrete  time  random  walk  X(t)  that  starts  at  a  given  point  Xq  and 
assume  that  over  a  given  time  interval,  At,  the  value  of  X(t)  can  increase  by  some  amount,  Ax  ,  decrease 
by  Ax,  or  remain  the  same.  If  we  start  at  a  given  time  t  and  population  x,  the  next  time  step  is  given  by: 

(10)  X(t  +  At)  e  {X{t)  +  Ax,  X(t)  -  Ax,  X{t)} 

We  refer  to  these  population  changes,  either  positive  or  negative,  as  spatial  changes,  and  each  of  these 
possibilities  of  spatial  change  has  a  probability  associated  with  it.  Define  P+  (■ t )  to  be  the  probability  of 
a  positive  spatial  change  over  a  time  step  At  starting  from  time  t,  define  P_(t)  to  be  the  probability  of  a 
negative  spatial  change  over  a  time  step  At  starting  at  time  t,  and  define  Po(t)  to  be  the  probability  of  no 
change  over  a  time  step  At  starting  at  time  t.  In  other  words,  if  we  let  P(e)  denote  the  probability  of  the 
event  e  then 

P+(t)  =  P(X(t  +  At)  =  X(t)  +  Ax) 

(11)  P_(t)  =  P(A(t  +  At)  =  X(t)  -  Ax) 

P0(t)  =  P(X(t  +  At)=X(t)) 

The  values  of  the  probabilities  are  chosen  such  that  they  match  the  differential  equation  of  the  birth  death 
model,  and  are  given  by  P+(t)  =  bX(t)At,  P-(t )  =  dX(t)At,  and  Po(t)  =  1  —  (P+(t)  +  P_(t)).  Note  that 
we  cannot  a  priori  guarantee  that  P+,  P_  and  P0  are  a  probability  distribution.  However,  we  proceed  with 
a  formal  calculation.  Let  us  assume  that  it  is  possible  to  choose  At,  in  such  a  way  that  P+  +  P_  is  less  than 
one,  and  as  such,  P+(t)  +  P_(t )  +  Po(t)  =  1.  These  values  were  chosen  so  that  the  expected  value  of  each 
step  satisfies  E[X(t  +  At)  —  X(t)]  =  (b  —  d) Ax.  We  will  see  this  later  on. 
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Since  this  model  is  discrete  and  the  value  of  X  at  each  time  step  only  depends  on  the  value  at  the 
previous  step,  we  know  that  if  X(t  +  At)  =  x,  there  are  only  three  possibilities  that  can  lead  to  this  state, 
namely: 

(1)  The  value  of  X(t)  increased  so  that 

X(t)=x  —  Ax  and  X(t  +  At)  =  x 

(2)  The  value  of  X(t)  decreased  so  that 

X(t)  =  x  +  Ax  and  X(t  +  At)=x 

(3)  The  value  of  X(t)  remained  the  same  so  that 

X(t)  =  x  and  X(t  +  At)  =  x. 

It  is  important  to  note  that  these  events  are  disjoint,  and  X(t  +  At)  =  x  is  the  union  of  these  events. 
Therefore,  the  probability  of  the  event  X(t  +  At)  =  x  is  equal  to  the  sum  of  the  probabilities  of  the  these 
three  events.  We  have: 


P (X(t  +  At)  =  x)  =  ¥(X(t  +  At)  =  x  fl  X[t)  =  x  —  Ax) 

+  P(X(t  +  At.)  =  x  D  X(t)  =  x  +  Ax) 
+  P(X(t  +  At)  =  xP)X(t)  =  x) 


(12) 

Since  P(A\B)  =  ,  it  follows  that  P(^4  Pi  B)  =  P(A\B)P(B),  and  because  of  this, 


(13) 


P (X(t  +  At)  =  x  fl  X(t)  —  x  —  Ax)  = 

P  (X(t  +  At)  =  x\X(t)  =  x  —  Ax)P(X(£)  =  x  —  Ax) 


From  above,  we  know  that  P(A(£  +  At)  =  x\X(t)  =  x  —  Ax)  is  the  same  as  P+(t),  and  similarly  for  P~(t) 
and  Po(t).  The  above  equations  then  reduce  to: 


(14) 


P  (X(t  +  At)  =  x)  =  P(X(t)  =  x  —  Ax)(P+(£)) 

+  P(A(£)  —  x  +  A  x)(P-(t)) 
+  P(X(t)=x)(P0(t)) 


Recall  the  definition  of  expected  value: 

Definition  3.1.1.  Consider  a  random  variable  X  that  has  possible  outcomes  xi, . . . ,  xn  with  correspond¬ 
ing  probabilities  pi, . . .  ,pn.  The  expected  value  of  X  is  given  by  E[X]  =  Y^i=iPixi- 

Define  Y(t)  :=  X(t  +  At)  —  X(t).  We  have, 


E[Y(t)\X(t)  =  x] 


Now  consider  the  infinitesimal  mean 


(x  +  Ax  —  x)bxAt  +  (x  —  Ax  —  x)dxAt. 

+  (x  —  x)(l  —  bxAt  —  dxAt) 

AxbxAt  —  AxdxAt  +  0(1  —  bx(At)  —  dx(At)) 
AxbxAt  —  AxdxAt 
( b  —  d)xAxAt 


lim 

— >-0 


E[Y(t)\X(t)  =x] 
At 


( b  —  d)x  Ax 
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We  can  do  the  same  computations  with  the  second  moment,  where  Y (t)  is  defined  as  above: 

E[Y(t)2\X(t)  =  x\  =  (x  +  Ax  —  x)2bX(t)At  +  (x  —  Ax  —  x)2dxAt 

+  (x  —  x)2(l  —  bxAt  —  dx)(At) 

=  ( Ax)2bxAt  +  ( Ax)2dx(At )  +  0(1  —  bxAt  —  dxAt) 

=  (6  +  d)x(Ax)2  At 

E[Y{t)2\X{t)  =  x\ 


At 


=  (b  +  d)  x(Ax)2 


Then  consider  the  infinitesimal  second  moment: 


..  E [Y{t)2X{t)=x]  n  „  ,2 

lim  L  '  - 1  =  (b  +  d)x(Ax)2 

At^o  At  '  y  ’ 

3.1.2.  A  Partial  Differential  Equation  associated  with  the  Random  Walk.  Define  v(x,t)  to  be 
P(A”(t)  =  a;|A(0)  =  Xo).  We  find 

v(x,  t  +  At)  =  v(x  —  Ax ,  t)P+(t )  +  v{x  +  Ax ,  t)P-(t)  +  v(x,  t)Po(t) 

Replacing  with  the  values  for  P+(t),  P-(t),  and  Po(t)  yields: 

v(x,  t  +  At)  =  v(x  —  Ax,  t)bX(t)At  +  v(x  +  Ax,  t)dX(t)At 
+v(x,t)(l  —  (b  +  d)X(t)At) 

Subtract  v(x,t)  from  both  sides,  divide  by  At,  and  recall  that  at  time  t,  X(t)  =  x: 

v(x,  t  +  At)  —  v(x,  t) 

(15) 


At 


=  v(x  —  Ax,  t)bx  +  v(x  +  Aa:,  t)dx 


+v(x,t)(b  +  d)(x) 


Define  a  new  function,  w  given  by  w(x,t)  =  xv(x,t)  and  rewrite  (15): 
v(x,t  +  At)  —  v(x,t) 


At 


=  bw(x  —  Ax,  t)  +  dw(x  +  Ax,  t)  +  (b  +  d)w(x,  t) 

=  b(w(x  —  Ax,  t)  —  w(x,  t)) 

+  d{w{x  +  Ax,  t)  —  w{x,  t)) 


From  Taylor’s  Theorem,  we  know  in  that  general  for  any  continuous  function  with  continuous  derivatives, 
denoted  /,  we  have 


(16) 


f{x  -  Ax)  =  f{x)  -  A xfxix)  +  fxxix)  +  0((Ax)3) 


The  term  0((Ax)3)  indicates  that  the  expansion  continues,  but  the  value  of  the  remainder  depends  on  (Ax)3 
and  is  close  to  0  when  Ax  is  small,  as  it  is  in  this  case.  So  we  can  expand  as  above  to  find, 

(Ax)2 

w(x  ±  Ax)  —  u;(x,t)  =  ±(Ax)wj(x,t)  +  — - — wxx{x,t)  +  0((Ax)3). 


The  right  hand  side  of  ( 16 )  then  becomes 


b((-Ax)wx(x,t)  +  ^Y~wxxix,t))  +  {d)(iAx)wx{x,t)  +  ^—y—wxx{x,t))  +  0((Ax)3) 
which  reduces  to  (d  —  b)(Ax)wx{x,  t)  +  (b  +  d){^^-wxxix,  t)).  From  this  position,  we  are  left  with: 


(17) 


;(x,  t  +  At)  —  x(x,  t) 

At 


=  id  —  b)i Ax)wx{x,  t)  +  (b  +  d) 


(Aa 


wxxix,t)  +  0((  Ax)3). 


3.2  The  Birth  Death  Model 


22 


We  note  that  this  equation  looks  very  similar  to  a  discretized  version  of  the  partial  differential  equation 
(f8)  ^  {{d ~ b)xv{x ’  *}) +  ((5 + d)xv{x ’ 

In  fact,  taking  the  limit  as  At  and  Ax  tend  to  zero  in  a  specific  way,  one  can  recover  this  equation  exactly. 
Additional  details  can  be  found  in  1  . 

3.1.3.  Transitioning  from  a  PDE  to  a  Stochastic  Differential  Equation.  From  1  ,  we  get  that 
for  small  At  and  Ax,  the  probability  distributions  in  the  partial  differential  equation  above  are  approximately 
the  same  as  the  probability  distributions  for  solutions  to  the  discrete  stochastic  process.  Therefore,  the 
discrete  stochastic  process  that  models  this  ordinary  differential  equation  is: 

(19)  dX(t)  =  (b  —  d)X{t)dt  +  \/(b  +  d)X(t)dWt 

Later  on  in  the  paper,  the  terms  Wt  and  dWt  will  be  defined  and  explained  in  more  detail,  but  for  now,  we 
can  consider  them  to  be  the  added  “randomness”  to  the  system. 

From  1 16] ,  we  know  something  similar.  If  we  are  given  a  discrete  system,  and  desire  a  stochastic 
differential  equation  of  the  form: 

(20)  dX(t)  =  F(X)dt  +  G(X)dWt, 

then  the  F  term  consists  of  the  infinitesimal  mean  from  the  discrete  system  and  the  G  term  is  the  infinitesimal 
second  moment.  Therefore,  if  we  desire  a  stochastic  differential  equation  for  the  three  component  model, 
the  best  place  to  start  is  to  discritize  the  system. 

3.2.  Brownian  Motion 

3.2.1.  Stochastic  Processes.  In  order  to  more  realistically  describe  situations  seen  in  real  life,  some 
models  include  stochastic  processes.  In  order  to  define  a  stochastic  process  we  begin  with  the  definition  of 
a  probability  space,  then  move  on  to  define  a  random  variable,  and  finally  a  stochastic  process. 

Definition  3.2.1  (Probability  Space).  A  probability  space  is  a  triple  (D,  U,P)  consisting  of 

(1)  A  sample  space  D,  which  is  a  set  that  consists  of  all  possible  outcomes. 

(2)  A  a-algebra  U.  A  cr-algebra  is  a  collection  of  subsets  of  D  that  contains  the  empty  set  0,  D,  and  is 
closed  under  countable  unions  and  complements. 

(3)  A  probability  measure  P  :  U  — >  [0, 1]  such  that  if  U\,U2,---  £  U  are  mutually  disjoint,  then 
P(u~i£4)  =  £~iP(t4). 

This  definition  of  a  probability  space  allows  us  to  define  a  random  variable: 

Definition  3.2.2  (Random  Variable).  Given  a  probability  space  (Q,U,P),  a  random  variable  is  a 
function  X  :  fi  — x  M  such  that  for  all  r  £  R,  the  set  {w|A(w)  <  r}  £  U 

We  can  view  a  stochastic  process  as  a  random  variable  indexed  by  time. 

Definition  3.2.3  (Stochastic  Process).  Given  a  probability  space  (D,[7, P),  a  stochastic  process  is  a 
function  I  :  !1  x  1  ->  I  such  that  for  each  i  £  R,  X(t)  =  Xt  is  a  random  variable. 

3.2.2.  Sample  Path.  With  this  in  mind,  we  can  define  a  sample  path. 

Definition  3.2.4  (Sample  path).  Let  X  :  !1  x  R  -y  K  be  a  stochastic  process.  For  each  w  £  the 
function  t  K >  X(ui,t)  is  called  a  sample  path. 

Since  Xt  consists  of  a  set  of  random  variables  indexed  by  time,  the  sample  path  of  Xt  for  u  would  be 
the  collection  of  values  Xt(uj).  This  is  a  function  from  time  to  the  possible  outcomes  of  Xt.  We  can  now 
define  a  continuous  stochastic  process. 

Definition  3.2.5  (Continuous  Stochastic  Process).  A  continuous  stochastic  process  is  a  stochastic  pro¬ 
cess  such  that  for  every  u>  £  D,  Xu(t)  is  a  continuous  function  for  t  £  [a,  6]. 
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3.2.3.  Brownian  Motion.  The  most  basic  stochastic  process  is  Brownian  motion,  which  is  also  known 
as  the  Wiener  process.  Originally  used  to  describe  the  motion  of  particles  suspended  within  a  fluid,  Brownian 
motion  is  named  after  the  botanist  Robert  Brown.  In  the  late  nineteenth  century,  he  observed  that  pollen 
floating  in  water  appeared  to  move  about  in  a  random  manner.  When  he  replaced  the  pollen  with  inorganic 
material,  he  noticed  that  the  motion  persisted.  Upon  plotting  the  motion  over  time,  he  noticed  that  it 
appeared  to  have  no  tangent  anywhere,  and  regardless  of  the  magnification  or  scale,  the  motion  appeared 
random.  This  motion  was  later  described  by  Albert  Einstein  to  be  caused  by  molecular  level  forces  on  the 
particle.  These  impacts  on  the  aggregate  cancel  each  other  out,  but  at  any  particular  time,  it  is  highly 
likely  that  one  side  of  the  particle  is  experiencing  more  force  than  the  other.  This  causes  the  particle  to 
move  about.  Mathematically,  this  process  was  described  first  by  Norbert  Wiener,  after  whom  the  process 
is  named.  We  abbreviate  the  Wiener  process  W(t)  or  Wt  for  each  t  £  [0,  T].  This  process  possesses  three 
important  and  distinctive  properties: 

(1)  W{ 0)  =  0 

(2)  W(t)  is  normally  distributed  with  mean  0  and  variance  t. 

(3)  For  any  sequence  0  =  to  <  t\,<  •••  <,tn  =  T,  the  random  variables  W(ti)  —  W(to),  Wfc)  — 
W(£i), . . . ,  W(tn)  —  W(tn-i)  are  independent. 

We  can  consider  Brownian  motion  to  be  the  limit  of  a  symmetric  random  walk  as  both  the  spatial  size 
and  time  step  go  to  zero,  or  perhaps  consider  it  to  be  the  integral  of  noise  independent  of  frequency. 

It  is  important  to  note  that  for  each  t,  W(t)  is  a  random  variable,  with  expected  value  0.  However,  if  we 
were  to  graph  several  different  sample  paths  of  W(t),  we  would  note  that  they  have  entirely  different  graphs. 

Figure  [3T] traces  out  different  sample  paths  for  Brownian  Motion  from  time  t  =  0  to  t  =  1000.  We  note 
that  the  spread  or  variance  in  each  path  grows  with  time.  Each  solid  black  line  represents  one  standard 
deviation  away  from  the  mean.  If  we  were  to  graph  more  sample  paths,  and  observe  at  any  given  time  their 
value,  we  would  see  that  the  distribution  of  the  values  at  that  time  closely  resembled  a  normal  distribution. 
This  is  visible  in  the  histogram  drawn  from  the  end  values. 


3.3.  Stochastic  Integration 

A  stochastic  differential  equation  (SDE)  is  an  equation  of  the  form 

dX(t)  =  F(X(t),t)dt  +  G(X(t),t)dWu 

where  F,  G  are  functions  and  Wt  is  Brownian  motion.  In  order  to  give  proper  meaning  to  the  above  equation 
we  would  need  to  develop  a  significant  amount  of  theory  that  is  outside  the  scope  of  this  project  and  a 
typical  undergraduate  syllabus. 

Instead  we  will  settle  for  a  more  heuristic  approach  and  provide  references.  The  manuscript  of  Evans 
4  is  particularly  good. 

We  can  gain  some  insight  into  the  SDE  by  discretizing  it: 

Examining  the  discrete  case,  we  see  a  stochastic  differential  equation  takes  the  form 


(21)  X(t  +  h)  =X{t)  +  hdX  =  X{t)+F(X(t),t)h  +  G(X(t),t)(Wt+h  -  Wt) 


Recall  from  the  definition  of  Brownian  Motion  that  the  term  Wt+h  —  Wt  is  a  normally  distributed 
random  variable  with  mean  0  and  variance  h.  A  stochastic  differential  equation  is  typically  given  in  the  form 
dX  =  Fdt  +  GdW ,  where  X  is  the  unknown  stochastic  process,  F  is  a  deterministic,  or  drift  function,  and 
G  is  a  random,  or  diffusion  term.  Both  F  and  G  can  be  functions  of  X (t)  and  f,  and  exist  in  the  following 
function  spaces: 


L1(0,T)  is  the  space  of  stochastic  processes  F  such  that  E  /0T  |F|dW 
L2(0,  T )  is  the  space  of  stochastic  processes  G  such  that  E  /QT  G2dW 


<  +00 
<  +00 
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Figure  3.1.  Several  sample  paths  of  Brownian  Motion  with  the  curves  x  =  y/t  and  x  =  -Vi 
outlining  the  standard  deviation  (top),  and  a  histogram  of  the  distribution  of  Brownian 
Motion  increments  (bottom). 


As  noted  in  4  ,  these  spaces  must  include  the  additional  technical  condition  that  the  stochastic  processes 
be  progressively  measurable.  We  omit  details  related  to  this  point  and  instead  refer  the  reader  to  |4|  for 
further  details. 

When  dealing  with  stochastic  differential  equations,  we  will  have  a  time  interval,  [0,  T],  and  see  terms 
such  as  udX  =  Fdt  +  GdW'\  where  X ,  F  and  G  are  stochastic  processes,  with  X,  F  £  L1([0,T])  and 
G  £  L2([0,T]).  This  dX  term  has  no  intrinsic  meaning,  but  rather,  implies  that  given  two  values,  s  and 
re[0,T], 

Fdt+  f  GdW 

J  r 

In  order  to  understand  what  this  actually  means,  we  should  break  the  differential  down  into  its  component 
parts.  The  left  side  of  the  equation  is  simple.  It  simply  subtracts  the  value  of  the  stochastic  process  X  at 


(22) 


A'OO  -  X(r)  = 
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time  r  from  the  value  of  X  at  s.  Remember  that  both  of  these  evaluations  are  random  variables.  The  right 
side  of  the  equation  consists  of  two  parts:  fr  Fdt  and  Jr  GdW.  We  shall  examine  each  of  these  separately. 

Recall  in  ordinary  calculus  how  integrals  were  introduced.  After  choosing  a  large  number  of  steps  n 
and  partitioning  the  interval  [r,  s]  into  to  =  r,  =  r  +  §Fr->  £2  =  r  +  2^r' ,  ■  ■  ■ ,  tn  =  s,  the  integral  was 
approximated  by  Riemann  sums  so  that 


(23) 


/a  n~1 

fdttt^f  (tk) 
k=o 


s  —  r 


n 


Something  similar  can  be  done  for  a  stochastic  process  F(t)  integrated  with  respect  to  time.  We  discretize 
the  time  interval  into  small  parts  and  evaluate  the  process  at  each  point  of  the  partition.  The  resulting 
integral  appears  identical  to  that  above 


n—  1 


(24) 


Fdt  ss  ^2  F(tk)  ' 


k= 0 


but,  contrasting  with  an  ordinary  integral,  however,  the  result  of  this  integration  is  now  a  random  variable, 
not  a  real  number.  Similarly  for  fr  GdW,  we  can  approximate  the  integral  as  defined  by  Evans  4  : 


(25) 


/a  n~1 

GdW  «  ^2G(tk)(Wi 

fc= 0 


wtJ. 


Note  that  this  is  also  a  random  variable. 

Just  as  in  ordinary  calculus  not  all  functions  F,  G  can  be  integrated  as  above.  Certain  technical  as¬ 
sumptions  must  be  made  about  the  adaptability  of  the  stochastic  process  X  to  the  filtration  generated  by 
Brownian  motion.  These  ideas  are  beyond  the  scope  of  this  project. 


3.3.1.  Basic  properties.  There  are  two  fundamental  properties  of  the  Ito  integral  that  we  need. 

The  first  property  tells  us  that  if  we  have  any  function  in  L1([0,  T])  integrated  with  respect  to  Brownian 
Motion,  then  the  expected  value  of  the  integral  is  zero.  In  other  words, 


(26) 


=  0 


The  second  is  the  Ito  Isometry.  The  Ito  Isometry  is  an  important  relation  between  dW  and  dt,  and  it 
permits  us  to  use  a  heuristic  that  is  vital  in  our  coding  methods. 


Theorem  3.3.1  (Ito  Isometry). 


(27) 


(  rT  \2 

rT 

E 

U FdW) 

=  E 

- 1 

•+0 

CM 

fcn 

_ _ O 

_ 1 

The  proof  is  given  in  |4  . 
namely 


Additionally,  an  important  special  case  of  this  theorem  arises  when  F  =  1, 


(28) 


E 


=T-0=T 


This  leads  to  the  heuristic  dW  S3  \fdt  or  dt  «  dW2,  which  will  be  used  for  various  calculations. 
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3.3.2.  Stochastic  Calculus.  Based  on  previous  discussion,  there  are  several  properties  of  stochastic 
differentials  worth  exploring.  The  first  is  that  the  differentiation  operator  is  linear.  So  given  two  stochastic 
processes,  X  and  Y,  then  we  can  compute  the  differential  of  any  linear  combination  d(aX+bY)  =  a  dX+b  dY 
where  R.  This  is  identical  to  what  we  see  in  ordinary  calculus. 

Contrastingly,  the  product  rule  for  stochastic  calculus  differs  from  ordinary  calculus.  Given  two  stochastic 
processes,  X1  and  X2,  such  that  dXi  =  F^dt  +  G\dW  and  dX2  =  F2dt  +  G2dW,  we  find 

(29)  d(X1X2)  =  XxdX2  +  X2dXx  +  GiG2dt. 

3.4.  Ito’s  Lemma 

With  everything  above,  we  can  state  Ito’s  Lemma.  Ito’s  Lemma,  which  is  also  known  as  Ito’s  Formula, 
is  one  of  the  most  important  results  in  stochastic  calculus  because  it  allows  us  to  find  solutions  for  stochastic 
differential  equations,  by  providing  an  analogous  version  of  the  Chain  Rule.  We  will  use  the  Ito  Isometry 
to  develop  a  partial  differential  equation  to  solve  for  the  probability  distribution  of  the  stochastic  three 
component  model  at  any  time. 

The  formal  statement  of  Ito’s  Formula  follows: 

Theorem  3.4.1  (Ito’s  Lemma).  Suppose  that  X  is  a  stochastic  process  with  differential 

(30)  dX  =  Fdt  +  GdW 

with  F  £  L1([0,T])  and  G  £  L2([0,T]).  Assume  that  a  function  u  :  R.  x  [0, T]  — >  R  is  continuous  with 
continuous  first  order  partial  derivatives  and  a  continuous  second  order  partial  derivative  .  Define  the 
stochastic  process 

(31)  Y(t):=u(X(t),t). 

The  stochastic  differential  for  Y  is  then 

dY  =  ^dt+^dX  +  ^dt 

=  (M  +  ^F  +  ^)dt  +  ^GdW 

proof  for  Ito’s  Lemma  is  actually  remarkably  simple,  and  an  outline  of  the  steps  involved  follows: 
By  using  the  stochastic  product  rule  and  mathematical  induction,  prove  Ito’s  Formula  for  functions 
of  the  form  Y{t)  :=  u(X(t),t)  =  ( X(t))n . 

Use  the  fact  that  all  functions  are  limits  of  polynomials  to  prove  Ito’s  Formula  for  functions  that 
are  polynomials  of  X(t)  or  f,  that  is  a  function  that  is  of  the  form  P{X{t))  =  cn(X(t))n  + 
cn-i(X  (t))n  1  +  •  •  •  +  ci(X(t))  +  c0. 

Use  the  product  rule  to  prove  Ito’s  Formula  for  separable  functions  of  X(t)  and  t ,  that  is,  for  a 
function  u(X(t),t)  =  f{X(t))g(t). 

Use  the  above  facts  to  show  that  if  any  sequence  of  polynomials,  Yn  =  un(X(t),t)  converges  to  a 
function  Y  =  u(X(t),t)  (with  the  derivatives  ^L,  and  likewise  converging  to 
and  fyO),  then  the  differentials  dYn  converge  to  dY. 

3.4.1.  Multi-dimensional  Ito’s  Formula.  Extending  Ito’s  Lemma  to  a  multidimensional  state  re¬ 
quires  some  additional  assumptions.  Assume 

(1)  W (t)  =  [Wi(t), . . . ,  Wn{t)\  is  an  n-dimensional  Brownian  motion.  Each  Wf  can  be,  but  is  not 
necessarily  independent  from  every  other  W). 

(2)  G (t)  is  an  Mnxm  valued  stochastic  process,  meaning  for  each  input  t,  an  n  x  m  matrix  of  random 
variables  is  the  output. 

(3)  F(f)  is  an  R”  valued  stochastic  process.  For  each  input  (t),  a  vector  of  length  n  of  random  variables 
is  the  result. 

(4)  F  and  G  are  in  the  spaces  [L1([0,T])]"  and  [L2([0,T])]”xm  respectively.  This  means  that  for 
each  index  chosen,  the  individual  process  that  generates  the  value  for  that  index  is  in  the  one 
dimensional  L-space. 


The 

(1) 

(2) 

(3) 

(4) 
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Then  we  can  say  that  an  Rra  valued  stochastic  process  X(f)  =  . . .  ,Xn(t)\  over  the  time  interval 

[0,  T]  satisfying  the  property 


(33) 


X(S)  -  X(r) 


Fdt 


GdW 


with  s,r  G  [0, T],  F  G  [L1([0, T])]ra,  and  G  G  [L2([0,  T])]raxn  has  the  stochastic  differential 
(34)  dX  =  Fdt  +  GdW 


which  means  that 

m 

(35)  dXi  =  Fidt  +  J^GijdWj  i  =  l,...,n. 

o= i 

With  these  definitions  in  mind,  we  can  state  Ito’s  Lemma  for  n-dimensional  stochastic  processes. 


Theorem  3.4.2  (Ito’s  Lemma  in  n  dimensions).  Suppose  that  X  is  a  stochastic  process  with  differential 


(36) 


dX  =  Fdt  +  GdW 


where  Fg  [L2([0,  T]]  ”  and  Gg  [L2([0,  T])] "xm. 
with  continuous  first  order  partial  derivatives 
=  l,...,n).  Then 


.  Assume  that  a  function  u  :  R"  x  [0,  T }  — >  Rn  is  continuous 
(i  =  l,...,n)  and  continuous  second  order  partial  derivatives 


(37) 


d(u(X(i),t)) 


o  n  ^  n  n 

O  il,  v — >  OU  1  v — r  v — > 


Ox., 


i= 1  3=1 


d2u 

dxidxj 


m 

EG^G^ 

i=i 


The  proof  follows  a  very  similar  line  of  reasoning  as  the  one  dimensional  formula  and  will  not  be  discussed 
here.  Instead,  see  [12  for  further  details. 


3.5.  The  Fokker-Planck  Equation 

Solving  stochastic  differential  equations  can  be  difficult  at  times,  while  solving  ordinary  differential  and 
partial  differential  equations  is  a  field  of  study  with  a  much  broader  base.  Therefore,  if  we  can  reduce  a 
stochastic  differential  equation  to  a  partial  differential  equation,  we  may  be  able  to  gain  additional  insight 
into  the  behavior  of  the  SDE. 

The  concept  that  draws  the  connection  between  the  PDE  and  the  SDE  is  the  notion  of  the  probability 
density  of  a  stochastic  process.  As  in  the  discussion  of  Ito’s  Lemma,  if  we  are  dealing  with  a  function  of  time 
and  a  stochastic  process,  the  probability  density  will  also  be  a  function  of  time  and  the  values  the  stochastic 
process  takes  on.  Before  we  begin,  it  is  important  to  distinguish  between  a  stochastic  process  X(t)  and  its 
probability  density  p(x,t).  We  know  that  at  any  given  time  t  for  which  the  stochastic  process  X  is  defined, 
evaluating  X  at  t  results  in  a  random  variable.  The  probability  density,  p{x ,  t)  at  this  same  time  t  gives  the 
distribution  of  the  random  variable  X(t). 

Now,  if  we  consider  the  stochastic  differential  equation 


(38) 

dX  = 

Fdt  +  GdW 

and  the  function 

(39) 

Y(t) 

=  u(X(t),t) 

with  the  following  conditions: 

(40) 

du  du 
dt  >  dx  5 

&  €Cg°(R) 

F  GL^O.T] 

G  GL2[0,T] 
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where  Cg° 


is  the  space  of  continuous  functions  with  an  arbitrarily  large  number  of  continuous  derivatives 
that  tend  to  zero  at  infinity,  then  we  know  that  from  Ito’s  Lemma: 


(41) 


dY  =  ^dt  +  ^dX  +  iyG^dt 


=  (dirt+tF+\a^)dt  +  ^xGdW 

If  we  apply  our  understanding  of  the  stochastic  differential  to  these  terms  and  recall  that  the  expected 
value  of  a  process  integrated  with  respect  to  dW  is  0,  the  following  is  a  result  of  taking  the  expected  value 
of  both  sides  of  the  equation  for  Y : 

(42) 


dt 


=  E 


du  du  1  d2u  2 
dt  +  faF+  28Y2 


Defining  the  probability  distribution  p(x,  t)  of  the  stochastic  process  X  as  above  and  replacing  the  expected 
value  term  with  its  definition  results  in 


(43) 


d  f°°  .  .  f°°  ( du  du  1  d2u  , 

-]_j(x,i)u(x,t)dx  =  j_j(x,t)  (-  +  -F  +  -—G- 


dx 


Rearranging  the  equation  and  changing  some  notation  simplifies  it  slightly 

r  ( d(pu 


(44)  7-00  V  dt 

Integration  by  parts  yields  the  following 

du 

pm 

J» 

/> ■£* 


-pIt )dx  =  0 


dt 


dx 


(45) 


d  (pu)  dp 

dt  ~  u~dt 

UPF l-oo  -  J 

d(pG2) 

~u  dx 


dx2 


ud{PF±dx 

dX 


u^dx 

dX2 


Since  u  and  € 


-<oo  ( 
" 0  V 


further  reduce  the  above  to 
(46) 


when  evaluated  at  infinity,  they  are  zero.  With  these  facts  in  mind,  we  can 
dp 


dx  =  0 


d{PF)  1  d2  (pG2 
dt  ^  dX  ~  2  dX2 
As  the  above  equation  will  hold  true  for  every  u  €  Cq°(R),  this  implies  that 


(47) 


dp(x,t)  d(p(x,t)F(x,t)) 


1  d2(p(x,t)G2(x,t)) 

2 


dt  dx  2  dx2 

This  final  equation  is  known  as  the  forward  Kolmogorov  (or  Fokker-Planck)  equation  for  probability 
distributions  of  solutions  to  stochastic  differential  equations.  This  particular  eciuation  is  one-dimensional. 
More  generally,  we  have  the  nrulti-dimensional  equation  as  follows,  with  the  following  stipulations: 

(1)  W (t)  is  an  n-dimensional  Brownian  motion 

(2)  G (t)  is  an  M”xm  real  valued  stochastic  process 

(3)  F(t)  is  an  M”  valued  stochastic  process 

(4)  F  and  G  are  in  the  spaces  [L1([0,T])]”  and  [L2([0,  T])]nxm  respectively. 

(5)  p(x,t)  is  the  probability  distribution  at  time  t  for  the  stochastic  differential  dX  =  F dt  +  GdWt. 
Then  p  satisfies  the  following  partial  differential  equation: 


(48) 


dp{x,  t) 
dt 


n  Gj  1  n  n 

=  t)pi  (*»  + 


i=  1 


2  '  dx2 

<=  1  3=1 


£)  ^  ^  GikG jk 


k= 1 
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With  sufficient  computing  power,  and  a  well-structured  numerical  method,  this  equation  would  allow  us 
to  supersede  the  Monte  Carlo  simulations  for  determining  the  end  probability  density  of  three  component 
SDE.  However,  since  it  would  require  large  amounts  of  computing  power  and  another  programming  language, 
this  particular  method  of  solution  was  bypassed. 
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CHAPTER  4 

Stochastic  Modeling 


4.1.  Creating  the  Stochastic  Model 

In  the  previously  described  birth  death  model,  a  differential  equation  was  derived.  From  this  differential 
equation,  we  were  able  to  create  a  stochastic  model.  This  was  done  by  noting  that  when  a  discrete  walk  of 
the  population  is  known,  the  infinitesimal  expected  value  and  second  moment  of  the  change  of  state  matrix 
provide  the  drift  and  variance  terms  included  in  the  stochastic  differential  equation.  In  this  section,  we 
present  the  analogous  derivation  for  (3CM). 


4.1.1.  State  Changes  and  Probability.  Recall  the  three  component  model: 


(3CM) 


dT 
dt 
dJ 
dt 
dV 
■  dt 


X  -  pT-  kTV 
kTV  -  SI 
pi  —  cV 


The  first  step  in  creating  a  stochastic  version  of  this  model  is  to  discretize  the  system.  This  is  done  using 
random  walks  on  the  state  space  of  values  of  T,  I,  and  V.  We  start  by  defining  three  discrete  random  walks 
that  correspond  to  the  populations:  T(t)  will  refer  to  a  random  walk  on  the  population  of  healthy  T-cells, 
I{t)  will  refer  to  a  random  walk  on  the  population  of  infected  T-cells,  and  V{t)  will  refer  to  a  random  walk 
on  the  virus  population.  Since  we  are  imposing  a  discrete-valued  model,  each  of  the  populations  must  take 
on  a  non-negative  integer  value,  and  all  changes  across  any  time  interval  At  will  be  an  integer.  If  we  shrink 
At  sufficiently,  then  we  can  assume  that  over  the  time  period  At,  only  one  state  change  occurs,  and  that 
state  change  will  alter  the  populations  by  +1,-1  or  0. 

It  is  rather  simple  to  list  the  possible  state  changes  of  the  vector  [T(t),  /(t),  V (t)]T  and  match  them  with 
their  respective  biological  meaning.  Each  state  change  will  be  denoted  by  A  Si  =  [a,  b,  c]T ,  where  a,  6,  and 
c  take  values  of  either  +1,  —1  or  0.  The  following  table  gives  all  the  possible  state  changes,  their  value,  and 
their  biological  interpretation: 


A  Si 

[a,  6,  c]T 

Biological  Interpretation 

A  Si 

[i,o,  op 

Production  of  a  Healthy  T-cell 

as2 

[-1,0, 0]T 

Death  of  a  Healthy  T-cell 

as3 

[-1,1,  of 

Infection  of  a  Healthy  T-cell 

A  S4 

[0,-1, 0]T 

Death  of  an  Infected  T-cell 

ASs 

[0,0, 1]T 

Production  of  a  Virus 

as6 

[0,0, -1]T 

Destruction  of  a  Virus 

AS0 

[0,0, 0]T 

No  Change 

In  a  manner  similar  to  that  of  the  birth-death  model  given  earlier,  each  of  these  possible  state  changes 
has  an  associated  probability,  derived  from  the  differential  equations  of  the  three  component  model.  They 
are  given  as  follows: 
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ASi 

Pi 

AS 1 

AA  t 

as2 

pT(t)At 

A  S3 

kT(t)V(t)At 

as4 

5I(t)At 

A  S5 

pI(t)At 

A  S6 

cI(t)At 

f) 

AS0 

i-Eft 

2=1 

4.1.2.  Expected  Value  and  Second  Moment.  It  is  known  that  the  expected  value  of  a  random 
variable  is  given  by  the  average  of  each  of  the  possible  values  the  variable  can  assume,  weighted  by  their 
corresponding  probabilities.  If  there  are  n  possible  events,  and  the  probability  of  event  Y,  is  pi ,  then 
E[Y]  =  In  this  case,  it  is  necessary  to  find  the  expected  value  of  the  change  of  state  vector,  and 

we  do  this  below.  Note  that  similar  to  the  calculation  in  the  birth-death  process,  there  are  parallels  between 
the  expected  value  and  the  ordinary  differential  equation.  If  we  divide  by  At  in  the  expected  value,  we  have 
a  term  that  looks  like  the  derivative  on  the  left,  while  the  right  sides  of  the  equations  are  the  same. 


(49) 


6 

E[AS|T(t)  =  a,  I(t)  =  (3,  V(t)  =  7]  =  J>*A  St 

i= 0 


A  —  pa  —  kaj 
kaj  —  5(3 
p(3  —  cq 


At 


The  second  moment  matrix,  which  measures  how  closely  each  variable  in  the  system  changes  with  respect 
to  each  other,  is  given  in  a  similar  way  to  the  expected  value: 


6 

(50)  E  [(AS)(AS)T]  =  Y,Pi(ASi)(ASi)T 

2=0 

As  calculated  from  the  values  given  above: 


(51)  E  [(AS)(AS)T|T(f)  =  a,I(t)  =  (3,  V(t)  =  7] 


A  +  pa  +  kaj 
—ka"/ 

0 


—ka  7 
ka  7  +  5(3 

0 


0 

0 

p(3  +  07 


At 


4.1.3.  Cholesky  Factorization.  The  second  moment  is  included  in  the  stochastic  model  under  a 
square  root,  so  it  is  necessary  to  compute  the  matrix  equal  to  (E  [(A5')(AS')r|T(t)  =a,I(t)  =  (3,V(t)  =7])5- 
In  this  case,  we  used  a  method  of  factorization  known  as  the  Cholesky  Factorization. 

Theorem  4.1.1  (Cholesky  Factorization).  A  real,  symmetric  matrix  A  is  positive  definite  if  and  only  if 
A  can  he  decomposed  into  A=LLT ,  where  L  is  a  unique,  real,  lower  triangidar  matrix  with  strictly  positive 
diagonal  entries. 


Computationally,  this  matrix  is  not  difficult  to  find.  We  know  that  L  is  of  the  form: 


(52) 


L  = 


h. 


l3,l 


0  0 

1-2,2  0 

h,2  h,3 


Therefore, 
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4.2.  A  SDE  based  on  the  Discrete  Model 

Before  we  can  model  a  stochastic  process  for  HIV,  we  must  first  rigorously  define  such  a  process  and 
prove  it  to  be  representative  of  this  biological  phenomena.  We  already  have  a  clear  idea  of  the  situation  we 
seek  to  model.  This  comes  from  the  discrete  model  shown  earlier;  therefore,  we  base  the  properties  of  the 
stochastic  model  on  the  properties  of  the  discrete  model. 

Recall  from  the  previous  section  that  at  time  t,  the  state  change  AS  possesses  the  following  two  prop¬ 
erties: 

A  —  n a  —  kaj 

(1)  E[AS'|T(t)  =  a,  I(t)  =  f3,  V(t)  =  7]  =  ka'y  —  S/3 

p(3  —  c'y 
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A  +  pa  +  ka'y  —ka'y  0 

—ka'y  ka  7  +  6/3  0 

0  0  p/3  +  07 

We  desire  to  create  a  stochastic  process  with  the  same  properties. 


(2)  E  [( AS)(AS)T\T(t )  =  a,  I(t)  =  (3,  V(t)  =  7]  = 


4.2.1.  Infinitesimal  Mean  and  Second  Moment.  From  16],  we  know  that  we  can  create  a  stochas¬ 
tic  differential  equation  for  the  three  component  model  from  the  discrete  model  by  using  the  infinitesimal 
means  as  the  coefficients  for  the  dt  or  “drift”  term  and  using  the  second  moment  for  the  dWt  or  “diffusion” 
terms. 

In  the  previous  section,  we  calculated  the  expected  value,  or  mean  and  the  second  moment  of  the  discrete 


model.  From  this  point,  based  on  the  results  of  Tuckwel  and  LeCorfec  16  ,  we  can  use  the  infinitesimal 


mean  and  second  moment  to  generate  a  stochastic  differential  equation.  These  infinitesimal  mean  and  second 
moment  are  computed  by  taking  the  standard  mean  and  second  moment  calculations,  than  dividing  by  At 
and  taking  the  limit  as  At  goes  to  zero.  In  this  way,  they  are  computed  over  an  “infinitesimally  small” 
period  of  time. 

Recall  that  the  mean  was: 


(53) 


E[A5|T(t)  =  a,  7(f)  =  0,  V(t)  =  7]  =  £  Pi  AS,  = 


i= 0 


A  —  pa  —  ka'y 
ka'y  —  5(3 
p/3  —  07 


At 


Dividing  by  At  and  taking  the  limit  as  At  goes  to  zero  gives  us  the  infinitesimal  mean: 


(54) 


lim  ms\nt)  =  a,m  =  (3,m  =  ,]  =  = 


At — ^0 


At 


1=0 


A  —  pa  —  ka'y 
ka'y  —  S/3 
p/3  —  cy 


Similarly  with  the  second  moment  matrix,  we  can  calculate  the  infinitesimal  second  moment  to  be: 


(55) 


lim 
At — >-0 


E  [(AS)(AS)T\T(t)  =  a,  J(t)  =  /3,  V(t)  =  7] 
At 


A  +  pa  +  ka  7 
—ka'y 
0 


—ka'y  0 

ka'y  +  6/3  0 

0  p/3  +  07 


4.2.2.  Piecing  Together  the  Stochastic  Differential  Equation.  Since  we  know  stochastic  differ¬ 
ential  equations  are  of  the  form 

dX{t)  =  F(X,  t)dt  +  G(X,  t)dWt 

Where  F  is  the  infinitesimal  expected  value  of  the  state  change  and  G  is  the  infinitesimal  second  moment, 
we  have  the  stochastic  differential  equation  for  the  three  component  model: 


T 

A  —  pT  —  kTV 

d 

I 

= 

kTV  -  SI 

dt 

V 

pI-cV 

'y/X+Jt/T+WV 


—kTV 

+A TJtT+kTV 

0 


0 

yj{kTV  +  SI)  - 

0 


(-kTV)2 

X+fiT+kTV 


0 


0 


dWt 


y/pl  +  CV 


This  also  allows  us  to  state  the  multidimensional  Fokker-Planck  Equation  for  the  stochastic  three  com¬ 
ponent  model.  Since  we  know  that  F  and  G  satisfy 


X  —  pT  —  kTV 

'A  +  pT  +  kTV 

-kTV 

0 

kTV  -  SI 

G  = 

-kTV 

kTV  +  SI 

0 

1 - 

1 

O 

1 _ 

0 

0 

pi  +  cV 

(56) 


F  = 
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we  find  from  the  first  half  of  the  equation 


(57) 


3 

=  (A  -  pT  -  SI  +  pi  -  cV) 

i= 1 


From  the  second,  if  we  make  the  following  substitutions  for  G 


'A  +  pT  +  kTV 

-kTV 

0 

a 

c 

O' 

-kTV 

kTV  +  51 

0 

= 

c 

b 

0 

0 

0 

pi  +  cV 

_0 

0 

d 

then  we  see 


3  3  3 

GikGjk 

2—1  J  =  1  fe  =  l 


(a  +  c)2  +  (6  +  c)2  +  d2 

(A  +  pT  +  kTV  -  kTV )2  +  (fcTU  +  5/  -  kTV )2  +  (pJ  +  cU)2 

(A  +  pT)2  +  {SI)2  +  {pi  +  cV)2 

A2  +  2A  pT  +  {pT)2  +  (<5/)2  +  (p/)2  +  2  {pcIV)  +  (cV)2 


With  all  of  this  in  mind,  we  can  write  the  multi  dimensional  Fokker-Planck  equation  for  the  stochastic  Three 
Component  Model: 

bOM)(A  -  pT  -61  +  pI-cV)\ 

+  l^r2  [p(x,t){ A2  +  2A pT+(pT)2  +  {SI)2  +  {pi)2  +  2 {pcIV)  +  { cV )2] 

4.3.  Existence  and  Uniqueness  of  Solutions  to  SDEs 


The  end  goal  of  crafting  the  SDE  is  to  find  the  solution  to  the  stochastic  differential  equation,  or  at 
the  very  least,  identify  its  important  properties.  Since  our  model  is  nonlinear,  it’s  quite  difficult  to  obtain 
an  analytic  representation  of  a  solution,  if  one  even  exists,  so  we  resort  to  computational  means  in  order 
to  calculate  sample  paths  of  solutions.  Regardless  of  how  we  attempt  to  obtain  a  solution,  it  is  important 
to  verify  that  the  solution  we  are  searching  for  exists  and  is  unique.  For  this,  we  can  use  theorems  from 
and  Mao  11. 


0ksendal 


12 


Theorem  4.3.1  (Mao).  Assume  that  the  stochastic  process  X{t)  satisfies 


(58) 


dX{t)  =  f{t,X(t))dt  +  g{t,  X(t))dWt 


with  the  given  initial  condition  A(0)  =  Xq  £  [L2(Md)]n;  where  f  and  g  satisfy  the  following  conditions  for 
some  T  >  0: 

(1)  (Local  Lipschitz)  There  exists  C  >  0  and  a  compact  K  C  such  that  for  all  t  £  [0,  T]  and  all 
x,y  £  K 

max.  {\f  (fix)  -  f(t,y)\2,  \g(fix)  -  g(t,y)\2}  <  C\x  -  y\2 

(2)  (Stochastic  Monotonicity)  There  exists  C  >  0  such  that  for  all  (t,x)  £  [0,  T\  x  Rd  such  that 


xTf(t,x)  +  -\g(fix)\2  <  C{l  +  \x\2). 


Then,  there  exists  a  unique  solution  X(t)  to  (58)  in  [L2([0,T];1 


For  the  SDE  that  we  have  derived  and  stated  at  the  end  of  the  previous  section,  the  quadratic  nature  of 
the  first  and  second  moments  and  bounds  from  the  previous  chapter  allow  /  and  g  to  satisfy  both  the  local 
Lipschitz  condition  and  the  Stochastic  Monotonicity  condition.  Hence,  we  are  guaranteed  the  existence  of  a 
unique  solution  on  some  time  interval,  though  we  do  not  control  the  extent  of  time  over  which  the  solution 
is  known  to  exist.  A  more  difficult  question  that  remains  unanswered  is  the  global  existence  of  this  solution. 
That  is,  can  this  solution  be  continued  to  an  arbitrarily  large  time  interval,  say  [0,  T]  for  any  T  >  0,  or  is 
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their  a  maximal  time  of  existence?  This  question  is,  unfortunately,  far  outside  the  scope  of  the  thesis,  and 
hence  we  will  not  pursue  it  further. 

Of  course,  from  a  modeling  or  applications  perspective,  explicitly  determining  the  solution  would  be 
much  more  helpful  and  provide  an  enormous  amount  of  information.  Since  tools  for  deriving  such  solutions 
for  SDEs  of  our  type  are  not  known,  we  proceed  instead  by  computing  solutions  numerically.  The  next 
chapter,  then,  focuses  on  numerical  methods  for  approximating  solutions  to  stochastic  differential  equations, 
and  more  specifically,  our  three-dimensional  stochastic  model. 
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CHAPTER  5 

Numerical  Methods 


5.1.  Numerical  Methods  of  Computation 

Solutions  to  stochastic  differential  equations  are  difficult,  if  not  mathematically  impossible  to  determine. 
Thus,  extended  amounts  of  computation  must  be  done  in  order  to  calculate  the  solution.  Methods  exist  to 
approximate  solutions  to  these  equations,  both  over  intervals  of  time  and  at  specific  points  of  time.  One 
of  the  most  basic  ones  is  Euler’s  Method.  This  particular  method  is  part  of  a  larger  group  of  methods  of 
approximations  known  as  finite  difference  methods. 

5.1.1.  Euler’s  Method.  Euler’s  method  is  a  way  of  solving  differential  equations  with  specific  initial 
values.  Consider  a  function,  x  over  the  time  interval  [to,tf],  with  its  derivative,  x' ,  expressed  as  a  function 
of  x  and  time:  x'  =  f(t,x(t))  and  an  initial  value,  x(t o)  =  Xq.  Then  divide  up  the  time  interval  into  n  steps 
of  size  h ,  such  that  tk  =  to  +  kh  and  tf  =  to  +  nh  Then  for  each  step,  x{tk+ i)  can  be  approximated  by 
x(t.k)  +  h  ■  f{tklx(tk)). 

At  the  end  of  this  these  calculations,  the  error  between  the  actual  value  of  the  function  and  the  computed 
value  of  the  function  is  proportional  to  the  step  size,  h.  Therefore,  decreasing  h  will  decrease  the  error  by 
a  proportional  amount.  However,  doing  so  will  increase  the  number  of  computations  necessary  to  complete 
the  method  over  the  entire  time  interval.  Considering  the  time  for  each  calculation  is  typically  on  the  order 
of  seconds,  a  slight  increase  in  time  required  might  not  seem  like  a  significant  issue.  However,  when  the  same 
calculation  is  required  for  hundreds  or  thousands  of  trials,  the  impact  grows  that  much  larger. 


5.1.2.  Discretization  of  the  System.  Suppose  we  have  a  system  of  n  stochastic  differential  equations 
given  by 

dX  =  [dXu . . . ,  dXn]T  =  [h{t,  Xi(t), . . . ,  Xn{t)), . . . ,  fn{t,  Xi(t), . . . ,  Xn(t))f 


with  the  initial  conditions 

X(0)  =  [X1(0),...,X„(0)]T, 


and  we  desire  to  model  the  system  and  acquire  a  solution  for  the  time  interval  [0,T].  The  first  step  would 
be  to  discretize  time.  Suppose  our  desire  for  accuracy  of  the  model  within  the  constraints  of  the  system 
lead  us  to  require  k  time  steps.  We  are  then  left  with  the  sequence  of  time  values  to  =  0,ti  =  = 

tP _ ,tk  =  T.  A  commonly  used  notation  is  to  define  At  :=  so  the  sequence  of  time  values  becomes 

f0  =  0,  ti  =  At,  t2  =  2At. . . . ,  tk  =  kAt. 

From  here,  it  is  a  simple  matter  to  calculate  the  vector 


dX  o  =  [/i(0,  Xl(0),...,  X„(0)), . . . ,  /„(  0,  Xl(0),...,  X„(0))f 


The  approximation  for  X(fi)  is  given  by  X(ti)  =  X(fo)  +  At.  ■  dX o.  With  X(H)  now  known,  we  repeat  the 
process  for 

dX  i  =  . .  .,xn(ti)), . . .  ,fn(ti,xi{ti), . . .  ,x„(ti))]T 

and  use  it  to  calculate  X(f2)  =  X(ti)  +  At  ■  dX i.  This  process  is  repeated  through  X(tk)  and  the  end  result 
is  a  sequence  of  vectors  [X(fo),  •  ■  • ,  X(ffc)]  corresponding  to  the  given  sequence  of  time  values.  These  values 
can  be  plotted  to  give  an  idea  of  what  the  solution  looks  like,  or  used  computationally  for  other  purposes. 
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5.1.3.  Application  to  SDEs.  As  seen,  when  we  desire  to  model  the  solution  to  an  ordinary  differential 
equation,  the  computations  are  quite  simple  and  deterministic.  Once  you  go  through  them  once,  you  can  be 
assured  that  the  same  results  will  occur  every  time.  With  a  stochastic  differential  equation,  the  fundamentals 
are  quite  the  same.  The  time  interval  and  function  are  discretized,  and  the  formula  for  calculating  the  next 
value  of  the  equation  from  the  previous  one  is  exactly  the  same.  However,  when  the  discretized  derivative 
for  a  stochastic  process  is  calculated,  it  involves  a  “dW”  term.  In  the  discritized  sense,  we  recall  that 
dW  ss  W(tk+i)  —  W(tk)  and  that  each  difference  in  W  is  normally  distributed  with  mean  0  and  standard 
deviation  equal  to  the  time  difference.  Therefore,  we  know  dW  ~  W(tk+ i)  —  W(tk)  ~  N(0,tk+ 1  —  tk)-  Since 
most  computational  software  includes  a  pseudo-random  number  generator,  we  can  calculate  a  dW  whenever 
its  necessary.  The  end  result  of  one  run  of  Euler’s  method  will  then  be  one  sample  path  of  the  stochastic 
process  that  solves  the  stochastic  differential  equation. 

The  result  of  solving  a  stochastic  differential  equation  is  a  stochastic  process,  which  yields  a  random 
variable  when  evaluated  at  a  given  point.  Therefore,  if  we  were  to  run  two  separate  trials  of  Euler’s  method 
for  the  same  SDE,  we  are  very  likely  to  get  two  entirely  different  sample  paths.  It  is  necessary  to  run  many 
such  trials  in  order  to  determine  the  probability  distribution  of  the  process  over  the  entire  interval  and  at  a 
specific  point. 


5.1.4.  Example  SDE.  The  following  figure  (|5.1[>  is  an  example  of  Euler’s  method  applied  to  the  process 
Y  =  W(t),  with  the  corresponding  SDE  dY  =  dW  and  initial  conditions  Y(0)  =  0.  We  will  consider  this  over 
the  interval  [0, 1]  with  step  size  .001.  We  will  use  the  properties  dW  ~  N( 0, 0.001)  to  aid  in  our  calculations. 
The  approximation  from  Euler’s  method  will  then  be  Y(tk+ i)  =  Y(tk )  +  At  •  dW . 


Figure  5.1.  Plots  of  Euler’s  Method  with  10  trials 

As  we  can  see,  the  sample  paths  begin  clustered  very  closely  around  the  line  y  =  t,  but  as  time  increases 
they  grow  more  and  more  spread  out,  and  resemble  a  normal  distribution  at  the  end.  (Figure |5.2[) 
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Figure  5.2.  Histogram  of  end  results  with  1000  trials 

The  time  required  for  this  rather  simple  stochastic  process  is  much  greater  than  for  the  ordinary  Euler’s 
method.  Generation  of  a  random  variable  is  a  relatively  intensive  computational  task,  and  requires  non-trivial 
amounts  of  time. 


5.2.  Boundary  Conditions 

When  dealing  with  our  particular  model  for  HIV  population  growth  within  the  human  body,  we  must 
make  some  choices  regarding  our  boundary  conditions,  in  order  to  keep  our  model  in  a  reasonable  realm  of 
biological  accuracy. 


5.2.1.  Definitions.  Boundary  conditions  come  into  play  when  the  process  being  modeled  nears  a  bound 
of  the  values  the  process  can  take  on.  In  our  model,  these  boundaries  are  drawn  from  the  biological  interpre¬ 
tation  of  the  model,  specifically,  they  state  that  none  of  the  populations  can  drop  below  zero.  This  makes 
sense,  since  it  is  impossible  to  have  negative  amounts  of  a  population.  Therefore,  since  the  computation  of 
the  model  does  permit  the  populations  to  drop  below  zero,  we  must  find  some  way  of  correcting  this.  In  the 
code,  these  are  triggered  with  if-statements  that  come  into  play  when  a  population  goes  below  zero. 

Once  the  population  is  negative,  there  are  three  ways  of  dealing  with  it: 

(1)  Set  the  population  to  the  boundary  level.  These  boundary  conditions  are  known  as  absorbing,  since 
if  the  population  drops  below  the  boundary  value,  the  boundary  “absorbs”  it. 

(2)  Set  the  population  to  the  positive  value  with  the  same  magnitude  as  the  negative  value.  These  are 
reflecting  boundary  conditions,  since  it  is  as  if  the  population  bounces  off  the  boundary. 

(3)  Set  the  population  to  some  random  value  between  the  two  above  options.  This  is  a  stochastic 
boundary  condition,  since  it  randomly  chooses  something  between  the  other  two. 


5.2.2.  Boundary  Conditions  within  the  Model.  This  figure  (5.3)  demonstrates  what  these  bound¬ 
ary  conditions  would  look  like.  Suppose  we  have  two  paths,  red  and  blue,  and  we  desire  to  have  a  boundary 
at  y  =  0.  If  we  assign  to  the  red  path  a  reflecting  boundary,  when  the  path  hits  the  boundary,  as  opposed 
to  continuing  on  its  path  (the  dotted  line),  it  returns  from  the  direction  it  came.  At  this  point,  the  absolute 
value  of  the  red  solid  line  and  the  red  dotted  line  are  the  same.  On  the  other  hand,  the  blue  line  has  an 
absorbing  boundary  condition.  When  the  blue  path  impacts  the  boundary,  the  path  takes  the  boundary 
value  and  stays  at  that  point  for  the  remainder  of  the  time. 

Coding  boundary  conditions  within  the  model  is  not  complex,  but  the  behavior  of  the  system  can  be 
changed  dependent  on  the  type  of  boundary  used.  Before  delving  into  the  effects  the  various  boundary 
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conditions  have  on  the  system,  it  is  important  to  note  that  the  T-cell  population  in  the  model  is  never 
actually  allowed  to  hit  zero.  Rather,  the  value  eps  in  MATLAB  is  used.  This  variable  eps  is  short  for 
epsilon,  and  has  a  value  of  approximately  2.2204  x  10-16.  This  near  zero  value  is  used  in  place  of  zero 
because  using  zero  would,  in  some  cases,  result  in  abnormal  behavior  in  the  calculations  within  the  code  or 
cause  undesired  premature  equilibriums  in  the  model. 


Figure  5.3.  Boundary  Conditions 

Reflecting  (Red)  and  Absorbing  (Blue)  Boundary  Conditioi 


When  simulating  the  equations,  both  reflecting  and  the  modified  absorbing  boundary  conditions  were 
applied  in  the  way  described  above.  The  two  figures  (5.4|5.5)  demonstrate  the  results  of  each  of  these  runs. 
The  important  thing  to  note  in  this  case  is  the  bimodal  distribution  of  the  healthy  T-cells  in  the  absorbing 
boundary  condition  figure.  Some  of  the  population  spikes,  then  drops  off  rapidly,  while  another  significant 
amount  of  the  population  spikes  than  slowly  fades  away.  This  is  in  direct  contrast  to  the  reflecting  boundary 
conditions,  which  all  see  the  spike  in  population,  followed  by  the  rapid  drop  off. 
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Figure  5.4.  Plots  of  3CM  (Absorbing  Boundary  Conditions) 
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Figure  5.5.  Plots  of  3CM  (Reflecting  Boundary  Conditions) 


5.3.  Coding  the  Three  Component  Model 

The  simulation  of  the  Three  Component  Model  was  developed  in  MATLAB,  and  there  were  several 
different  paths  taken  to  get  final  probability  distributions,  and  each  of  them  differ  in  their  method  of  com¬ 
putation. 

5.3.1.  The  Discrete  Simulation.  This  code  simulates  the  discrete  model  given  in  the  previous  section. 
As  with  all  of  the  code,  the  basic  idea  is  to  determine  how  the  populations  change  over  a  time  interval,  and 
run  that  simulation  many  times  in  order  to  obtain  data  that  can  be  analyzed  in  a  Monte-Carlo  method.  We 
begin  here  by  defining  the  number  of  trials,  the  length  of  the  time  interval,  and  the  number  of  steps  across 
the  time  interval.  We  also  define  the  state  change  vectors.  Recall  from  before  that  these  vectors  take  the 
form  [x,y,  z],  where  x,  y ,  and  z  can  only  take  values  of  +1,  —1,  or  0. 

Once  these  beginning  steps  are  complete,  initial  conditions  are  generated  for  the  first  trial.  From  these 
populations,  the  probability  of  each  state  change  vector  occurring  is  calculated.  Recall  that  these  probabilities 
are  tied  to  the  differential  equation  and  depend  on  the  population  values.  Next,  a  uniform  random  number 
from  the  interval  [0,1]  is  generated  by  MATLAB,  and  scaled  to  the  interval  [0,P],  where  P  is  the  sum  of 
the  probabilities.  This  interval  is  broken  down  into  subintervals  of  the  form  [0,pi),  [pi,pi  +  P2),  ■  ■  ■ ,  \pi  + 
■  ■  ■  +  Ps,  P).  In  this  way,  the  width  of  each  interval  is  directly  proportional  to  a  probability.  We  determine 
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Figure  5.6.  Flow  chart  representation  of  the  simulation  process 


which  interval  the  random  number  fits  into  and  select  the  corresponding  state  change  vector  to  add  to  the 
population.  This  constitutes  one  time  step. 

Once  the  entire  time  interval  has  been  crossed,  the  path  is  plotted  and  the  end  values  stored.  The  process 
is  repeated  for  each  trial,  and  at  the  end  of  all  of  the  trials,  the  average  plot  is  created,  as  is  a  histogram  of 
all  the  end  values.  In  this  way,  we  can  determine  the  percentage  of  cases  in  which  the  end  virus  level  was 
beneath  detectable  concentrations. 

The  flow  chart  (figure  |5.6[)  gives  a  general  idea  of  what  is  being  accomplished  in  the  code.  Structurally, 
the  code  contains  two  nested  for— loops,  one  for  the  time  interval  and  one  for  the  set  of  trials. 

5.4.  Coding  Stochastically 

Modeling  the  stochastic  version  offer  advantages  not  seen  in  other  forms  of  modeling.  Primarily,  when 
we  model  stochastically,  we  get  to  observe  the  path  each  trial  takes,  as  opposed  to  just  observing  the  average 
path.  With  this,  we  get  a  much  better  idea  of  the  distribution  of  paths  at  a  given  sample  time.  The  downside, 
however,  is  that  this  modeling  is  computationally  expensive  and  requires  significantly  more  trials  per  sample 
to  get  acceptable  results. 

5.4.1.  The  Basic  Stochastic  Model.  The  basic  stochastic  simulation  uses  MATLAB  to  create  thou¬ 
sands  of  sample  paths  of  the  populations  over  time,  with  each  path  using  the  same  set  of  coefficients  and 
initial  conditions.  The  actual  code  bears  more  resemblance  to  the  discretized  simulation  than  the  deter¬ 
ministic  one,  due  to  the  method  of  implementing  the  simulation.  For  each  sample  path,  or  trial,  the  time 
interval  is  divided  in  separate  intervals.  At  the  endpoint  of  each  interval,  the  expected  change  in  population 
is  calculated.  This  expected  value  change  is  the  exact  same  as  would  have  been  seen  in  the  deterministic 
model.  The  stochastic  model  differs,  however,  with  the  addition  of  a  variance  term.  The  variance  of  the 
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population  change  is  calculated,  and  a  random  vector  is  generated.  The  variance  matrix  is  multiplied  by 
the  random  vector  and  added  to  the  expected  value.  This  new  population  change  is  then  added  to  the  old 
population.  Once  this  entire  process  has  been  completed  over  the  entire  time  interval,  it  is  repeated  for  a 
set  number  of  trials.  At  the  end  of  all  the  trials,  we  see  the  average  sample  path  and  the  histogram  for  the 
end  populations. 


Figure  5.7.  A  simulation  of  the  stochastic  model  with  R=8.18 


Again,  we  note  that  the  green  line  represents  the  virus  threshold  for  detection,  and  that  the  majority  of 
paths  and  the  average  sample  path  are  below  that  detection  threshold.  However,  it  is  incredibly  important 
to  note  that  not  all  sample  paths  fall  below  there. 
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CHAPTER  6 

Conclusion 


6.1.  Results 

To  conclude  the  thesis,  we  will  summarize  the  theoretical,  computational,  and  biological  results  obtained 
within.  Details  can  be  found  in  Chapters  2  —  5,  but  for  brevity  and  completeness,  we  highlight  them  below. 

6.1.1.  Theoretical  Results.  The  theorems  obtained  in  previous  chapters  provide  information  about 
the  model  from  a  mathematical  perspective.  In  addition  to  verifying  facts  necessary  to  complete  the  math¬ 
ematical  analysis  of  the  system,  these  theorems  demonstrate  that  the  model  is  biologically  relevant,  and 
remains  so  for  all  times,  and  yield  insight  into  necessary  methods  of  simulation.  Below  is  a  brief  summary 
of  the  main  theorems. 

Theorem  1  (Existence  and  Uniqueness  of  Solution  for  the  Three  Component  Model).  Given  initial 
conditions  T( 0),  1(0),  and  U(0)  for  the  Three  Component  Model,  there  exists  a  unique  solution  T(t),  I(t), 
and  V  ( t )  that  satisfies  the  system  of  differential  equations  for  all  t  >  0  and  the  given  initial  conditions. 

This  theorem  proves  that  a  solution  exists  for  the  Three  Component  Model  and  that  this  solution  is 
unique  for  every  unique  set  of  initial  conditions.  This  represents  a  fundamental  step  in  the  analysis  of  this 
system  of  differential  equations.  While  the  three  component  model  has  been  used  previously  in  other  studies, 
no  article  within  the  literature  proves  this  fundamental  fact. 

Theorem  2  (Positivity  of  Solution  for  the  Three  Component  Model).  Given  initial  conditions  T(0)  >  0, 
1(0)  >,  and  U(0)  >0  for  the  Three  Component  Model,  the  unique  solution  remains  bounded  for  every  t>  0 
and  possesses  the  property  that  T(t)  >  0,  I(t)  >  0,  and  V(t)  >0  for  all  t>  0. 

With  the  above  result,  we  confirm  the  biological  relevance  of  the  model.  This  theorem  ensures  that 
given  positive  initial  values  of  each  population,  the  system  will  evolve  so  as  to  preserve  this  property,  and 
population  values  will  remain  positive  for  all  later  times.  Since  negative  population  values  are  a  biological 
impossibility  -  meaning  that  the  number  of  T-cells  or  virions  within  the  body  should  never  become  negative 
-  the  confirmation  that  this  event  is  also  a  mathematical  impossibility  within  our  model  helps  verify  the 
relevance  of  our  derivation. 

Theorem  3.  Dependence  of  the  equilibrium  on  the  reproductive  constant  R  Define  the  reproductive 
constant,  R  by 

n  _  kXP 
cSp. 

The  infective  equilibrium  is  stable  if  and  only  if  R  >  1.  In  this  case,  any  solution  of  the  three- component 
model  possesses  the  property  lim^+oo  V(t)  >  0  and  the  viral  infection  is  persistent.  Conversely,  the  non- 
infective  equilibrium  is  stable  if  and  only  if  R  <  1.  In  this  case,  any  solution  of  the  three- component  model 
possesses  the  property  lim^+oo  V(t)  =  0,  and  the  viral  infection  becomes  extinct. 

Finally,  this  last  theorem  demonstrates  that  our  analysis  of  the  differential  equations  can  be  based 
solely  on  an  analysis  of  the  coefficients,  and  in  particular,  the  value  of  a  single  parameter  R.  While  this  is 
mathematically  true,  the  biological  applicability  of  this  statement  may  fail  since  it  can  only  be  determined 
that  the  virus  population  either  persists  or  becomes  extinct  as  t  — >  oo.  The  theorem  is  unable  to  tell  us  at 
what  rate  this  occurs,  or  what  the  values  of  the  viral  population  are  for  earlier  time  periods.  One  might 
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guess  that  eventual  extinction  should  generally  imply  that  the  viral  population  remains  small  even  for  small 
times.  However,  as  we  will  see,  later  numerical  results  indicate  that  this  is  not  the  case. 

6.1.2.  Numerical  Results.  Recall  that  this  project  simulated  the  three  component  model  under  many 
different  conditions.  The  results  were  discussed  previously  in  the  thesis,  but  are  listed  here  for  convenience 
and  overview. 


Model 

Coefficients 

Sample  Paths 

R  <  1 

V(60)  <  40  Vffops 

ODE 

Random 

10000 

211 

2140 

SDE 

Mean  Valued 

10000 

0 

1417 

SDE 

Random 

10000 

300 

2536 

Table  6.1.  Values  from  the  differing  simulations  of  ODE  and  SDE  models 


Notice  that  the  proportion  of  trials  within  our  simulation  that  resulted  in  R  <  1  is  much  smaller  -  by 
nearly  a  factor  of  ten  -  than  those  which  resulted  in  viral  populations  less  than  40  virions  per  /iL  approximately 
two  months  after  the  initial  infection.  Hence,  a  standard  mathematical  definition  of  extinction,  namely 
limt_>00V(t)  =  0  which  (by  Theorem  |3j)  is  equivalent  to  R  <  1,  is  insufficient  to  accurately  describe  the 
behavior  of  the  viral  population  on  the  timescales  of  biological  relevance,  specifically  during  the  time  period 
up  to  a  few  months  after  initial  infection.  Instead,  a  more  precise  determination  of  viral  extinction  can  be 
made  by  measuring  whether  the  virus  population  will  remain  below  the  current  detectable  threshold  of  40 
virions  per  /xL,  and  under  this  measure,  the  probability  of  extinction  is  much  larger.  As  can  be  seen  by  Table 


6.1  the  probability  of  extinction  by  the  former  measure  is  only  around  2%,  while  the  probability  jumps  by 


a  factor  of  ten  to  approximately  20  —  25%  under  the  latter  notion  that  we  have  proposed.  In  some  types  of 
transmission,  such  as  passage  of  the  disease  from  a  mother  to  an  unborn  child  or  through  needlesticks  ,  this 
percentage  is  much  closer  to  current  estimates  of  the  probability  of  extinction  after  initial  infection  than  the 
criteria  R  <  1. 


Virus  Particles 


Figure  6.1.  Histograms  of  end  state  values  from  ODE  w/  RV  coefficients 
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Virus  Particles 


Figure  6.2.  Histograms  of  end  state  values  from  SDE  w/  mean  coefficients 
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Figure  6.3.  Histograms  of  end  state  values  from  SDE  w/  RV  coefficients 


The  histograms  (6.1  6.2  6.3)  display  the  end  population  values  for  each  of  the  healthy  T-cells,  infected 


T-cells,  and  virions,  as  described  by  a  particular  model.  The  top  set  of  figures  is  for  the  ordinary  differential 
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equations  simulated  with  random  variable  coefficients,  the  middle  set  of  figures  is  for  the  stochastic  differ¬ 
ential  equation  simulated  with  the  mean  valued  coefficients  for  each  of  the  parameters,  and  the  final  set  of 
histograms  is  for  the  stochastic  equations  simulated  with  random  variable  coefficients.  There  is  line  drawn 
on  the  histograms  for  the  viral  populations.  All  entries  to  the  left  of  this  bar  represent  those  simulations 
that  resulted  in  the  number  of  virus  particles  being  less  than  40  virions  per  fiL.  This  is,  as  mentioned  before, 
the  biological  threshold  for  persistence. 


6.2.  Conclusion 


Throughout  this  research,  we  sought  to  answer  the  question:  “How  does  creating  a  stochastic  model 
alter  results  obtained  from  simulating  the  Three  Component  Model?”  The  idea  for  the  question  came  from 
Tuckwell  and  Shipman.  They  found  17  in  their  simulations  of  the  model  of  ordinary  differential  equation 
that  viral  persistence  happened  at  a  rate  of  approximately  2  %.  This  was  done  by  calculating  randomly 
generated  R  values  and  counting  the  number  above  and  below  one. 

In  our  work,  we  did  three  things.  First,  we  redefined  the  notion  of  persistence.  We  note  that  the  end 
viral  population,  V(60),  could  less  than  the  threshold  of  detectability  permitted  by  modern  science  at  a 
much  greater  rate  than  expected  by  the  R  value.  What  does  this  mean?  It  indicates  that  the  methods  used 
by  Tuckwell  and  Shipman  are  not  as  accurate  as  could  be  hoped  for.  It  makes  clear  that  a  full  simulation  of 
the  model  is  required  to  get  accurate  data.  Therefore,  with  our  second  notion  of  extinction  and  persistence 
stemming  from  a  biological  perspective,  we  can  carry  out  the  simulations  again. 

With  full  simulation  of  the  ordinary  differential  equations,  we  discovered  that  viral  persistence  happened 
at  a  rate  of  roughly  ten  times  that  suggested  in  17  .  This  proves  rather  conclusively  that  biological  persis¬ 
tence  and  mathematical  persistence,  while  related,  are  not  identical,  and  simulation  is  necessary  to  realize 
the  full  implications  of  the  model.  In  reality,  persistence  can  happen  at  a  rate  of  up  to  90  %,  if  the  virus 
is  transmitted  by  blood,  or  around  1  %  via  intercourse  |5  .  While  2  %  from  17  seems  to  be  accurate, 


recall  that  these  coefficients  were  drawn  from  people  with  HIV  already,  and  therefore  are  biased  towards 
transmission. 

Second,  we  took  the  mean  valued  coefficients  and  simulated  the  stochastic  differential  equations  for  the 
Three  Component  Model.  These  coefficients  give  us  the  “average”  person,  in  a  sense,  and  the  stochasticity 
allowed  us  to  simulate  this  person  many  times,  with  slight  random  variations  each  time.  Based  on  the 
coefficients  alone,  we  would  expect  the  virus  to  persist  all  of  time,  but  as  we  demonstrated,  a  full  simulation 
would  be  required.  In  addition  to  this,  the  added  stochasticity  would  incorporate  the  randomness  inherent 
in  the  external  world.  The  virus  persisted  most  of  the  time,  but  far  from  all  of  it. 

Finally,  by  combining  the  first  and  second  methods,  we  simulated  the  stochastic  model  with  different 
sets  of  coefficients.  The  results  we  obtained  were  similar  to  those  of  the  first  method.  Persistence  happened 
nearly  eight  times  more  than  would  be  expected  from  the  coefficients  alone.  This  demonstrates  that  the 
stochasticity  also  alters  the  results  of  the  simulation  towards  a  more  realistic  end. 

Stochastic  models  tell  us  more  than  can  be  told  with  models  based  on  ordinary  differential  equations. 
We  learned  that  incorporating  stochasticity  noticeably  alters  the  results  we  would  otherwise  expect  from 
simulation.  This,  coupled  with  a  different  notion  of  persistence  shows  us  we  can  improve  the  model  to  more 
accurately  describe  real  world  conditions. 


6.3.  Areas  of  Further  Study 

Further  study  could  be  done  in  a  number  of  areas.  The  model  could  be  expanded  to  include  anti¬ 
retroviral  treatment  (HAART),  or  to  include  such  biological  phenomena  as  the  latent  cell  reservoir.  These 
models  would  then  more  closely  resemble  the  interaction  between  the  virus  and  the  body.  From  a  more 
mathematical  end,  the  Fokker-Planck  equation  described  previously  with  the  paper  could  be  solved  for  the 
probability  density  at  the  end  of  the  model.  Additionally,  the  stochastic  differential  equations  themselves 
could  be  studied,  and  analyzed  for  their  stability  and  equilibrium.  This  study  was  only  able  to  examine 
existence  and  uniqueness  of  solution. 
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While  the  first  set  of  proposed  extensions  describe  ways  to  expand  an  accurate  model  of  extinction 
probabilities  and  the  second,  methods  of  expanding  purely  mathematical  knowledge,  both  would  be  useful 
in  further  efforts  to  model  HIV  in  its  early  stages. 
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APPENDIX  A 

Code 


A. 0.1.  Deterministic  Code. 

ODEs  at  different  coefficient  sets: 


1  function  sol=PopDer (t, P) 

2  %P  represents  a  vector  of  initial  populations  in  the  form  [T,  I,  V] 

3 

4  lambda  —  4.3e— 2; 

5  mu  =  .020; 

6  k  =  .  19e— 3; 

7  A  =  .080; 

8  p=  98; 

9  c=  3 . 8 1 ; 

10 

ll  sol  =  [lambda— mu*P (1)—  k*P (1) *P (3) ; k*P (1) *P (3)— A*P (2) ;p*P (2)— c*P (3) ] ; 


1  function  sol=PopDer2 (t , P ) 

2  %P  represents  a  vector  of  initial  populations  in  the  form  [T,  I,  V] 

3 

4  lambda  =  .2; 

5  mu  =  .0043; 

6  k  =  4 . 8e— 3; 

7  A  =  .13; 

8  p=  7100; 

9  c—  2.06; 

10 

ll  sol  =  [lambda— mu*P (1)— k*P (1) *P (3) ; k*P (1) *P (3)— A*P (2) ;p*P (2)— c*P (3) ] ; 
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Generate  figures  and  statistics  from  deterministic  model: 


1 

clear 

2 

clc 

3 

cl  f 

5 

lambda  =  .11; 

6 

mu  =  .01089; 

7 

k  =  .001179; 

8 

A  =  .3660; 

9 

p=1426 .8; 

10 

c=3 .074; 

11 

12 

R= (k*lambda*p) / (mu*A*c) 

13 

14 

[t,  y]  =ode45  (SPopDer,  [0,60],  [200, 0,  4e— 7] ) ; 

15 

subplot  (2, 3, 1) ,  plot (t , y ( : , 1 ) ) ,  xlabel('Time 

(days )'), ylabel (' Healthy  T— cells  (units?)'); 

16 

subplot (2 ,  3, 2 ) ,  plot (t, y ( : , 2) ) ,  xlabel('Time 

(days) '), ylabel (' Infected  T— cells  (units?)'); 

17 

subplot  (2, 3, 3) ,  plot  (t, y ( : , 3) ) ,  xlabel('Time 

(days )’) , ylabel (' Virus  ... 

(part  ides /microliter)  '),  line  ([0,60],  [40 

, 40] , ' color ' , ' r ' ) ; 

18 

[t ,  y  ]  =ode4  5  ( @PopDer2 ,  [0,60],  [200,  0,  4e-7] ) ; 

19 

subplot (2 , 3, 4 ) ,  plot (t, y  ( : , 1) ) ,  xlabel('Time 

(days )'), ylabel (' Healthy  T— cells  (units?)'); 

20 

subplot (2, 3, 5) ,  plot (t , y ( : , 2 ) ) ,  xlabel('Time 

(days )'), ylabel (' Inf ected  T— cells  (units?)'); 

21 

subplot  [2, 3, 6) ,  plot  (t, y ( : , 3) ) ,  xlabel(’Time 

(days) ' ) , ylabel ( 'Virus  . . . 

( part icles/micro liter) '),  line ([0,60], [40 

, 40] , ' color ' , ' r ' ) ; 
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A. 0.2.  Deterministic  Code  with  RV  coefficients. 

Generate  truncated  normal  distributions 


1  function  z  =  truncatednormal (N, mu, sigma, trunc.low, trunc_high) 

2  %  truncatednormal:  trncated  normal  random  variable  generator 

3  % 

4  %  input  explanations 

5  %  N  —  size  of  resultant  matrix/array 

6  %  mu  —  mean  of  untrncated  distribution 

7  %  sigma  —  std  dev  of  untruncated  distribution 

8  %  trunc.low  —  low  truncation  point 

9  %  trunc_high  —  high  truncated  point 
10  % 

n  %  output  explanation 

12  %  z  —  array  of  size  N  of  truncated  random  variables  with  above  properties 

13  % 

14 

15  %  set  the  defaults 

16 

17 

18  %  If  N  is  not  entered,  set  the  default  value  of  N  to  be  1 

19  if  (nargin  <1)  |isempty(N) 

20  N=1 ; 

21  end 

22 

23  %If  no  mean  is  given,  set  it  to  be  0 

24  if  (nargin  <  2) |isempty(mu) 

25  mu=0; 

26  end 

27 

28  %If  no  standard  deviation  is  given,  set  it  to  be  1 

29  if  (nargin  <  3) | isempty (sigma) 

30  sigma=l; 

31  end 

32 

33  %If  no  low  end  of  truncation  is  given,  set  it  to  be  negative  infinity 

34  %  prob_low  is  the  probability  that  the  the  value  of  the  random  variable  is 

35  %  less  than  the  truncation  point  for  the  given  distribution,  if  trunc.low  is 

36  %  empty,  set  prob_low  to  be  zero 

37  if  (nargin  <  4) | isempty (trunc.low) 

38  trunc_low=—  inf ; 

39  prob_low=0; 

40  else 

41  prob_low=normcdf ( (trunc.low— mu) /sigma) ; 

42  end 

43 

44  %If  no  upper  end  of  trunation  is  given,  set  it  to  be  infinity 

45  %  prob.high  is  the  probability  that  the  value  of  the  random  variable  is 

46  %  less  than  the  truncaion  point.  If  trunc.high  is  empty,  set  it  to  be  one. 

47  if  (nargin  <  5) | isempty (trunc.high) 

48  trunc_high=inf ; 

49  prob_high=l; 

50  else 

51  prob_high=normcdf  (  (trunc.high— mu)  /sigma)  ; 

52  end 

53 

54  %  Test  for  correct  order  of  truncation  points 

55  if  t runc_low>trunc_high 
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56  error  'Must  have  trunc.low  less  than  or  equal  to  trunc_high' 

57  end 

58 

59  %Generate  a  size  N  standard  normal  random  variable 

60  r=rand(N); 

61 

62  %scale  to  the  interval  [plo,phi] 

63  r=prob_low+  (prob_high— prob_low)  *r; 

64 

65  %Invert  using  the  standard  normal 

66  z=norminv ( r ) ; 

67 

68  %Rescale  and  shift  using  given  mean  and  standard  deviation 

69  z=mu+sigma* z ; 
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Generate  coefficient  sets 


1  function  [lambda, mu, k, A, p, c]  =  CoefCreator; 

2  %  This  function  generates  random  coefficents  for  use  in  3CM. 

3 

4  %  We  create  the  random  coefficients  based  on  the  data  from  Stafford,  et  all 

5  %  The  function  "truncatednormal (N, mu, sigma, trunc.low, trunc_high) "  generates 

6  %  an  array  of  random  variables  of  size  N,  with  mean  mu  and  standard 

7  %  deviation  sigma,  truncated  to  be  in  the  interval  [ trunc.low,  trunc.high]  . 

8  %  The  interval  of  truncation  is  given  to  be  [0,inf) . 

9 

10 

11  %  Here  we  use  the  max  and  min  values  from  Stafford  et  all  as  the  truncation 

12  %  points. 

13 

14 

15  lambda  =  truncatednormal (1,  . 1089,  . 05727,  . 043,  .  20) ; 

16  mu  =  truncatednormal (1,  . 01089,  . 005727,  .  0043,  .  020) ; 

17  k  =  truncatednormal (1, . 001179, . 001422, . 00019, . 00480) ; 

18  A  =  truncatednormal (1,  .3660,  .193,  .13,  .80) ; 

19  p  =  truncatednormal (1, 1426 . 8, 2049 . 36, 98, 2697) ; 

20  c  =  truncatednormal (1, 3 . 074, 0 . 64, 2 . 06, 3 . 81) ; 

21  %c=3 ; 
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Ordinary  Differential  Equations: 


1 

2 

3 

4 

5 

6 

7 

8 


function  sol=Det _3CM (t , P , lambda, mu, k, A, p, c) ; 

%  This  function  is  the  ODE  for  the  deterministic  3CM. 

%  t  is  time,  P  is  population,  coef  is  the  coefficient  vector 


sol  =  [ lambda— mu* P (l)-k*P (1) *P (3) ;k*P (1) *P (3)-A*P (2) ;p*P (2)-c*P (3) ]  ; 
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Simulate  ODE: 

1 

function  Data=ODE_Simulator ( ) 

3 

clear 

4 

clc 

5 

6 

7 

tic 

%  Determine  the  number  of  trials 

8 

trials=100; 

10 

%  Create  a  holding  matrix  for  the  end  populations  and  the  R  value  [T 

I  V  R] 

11 

Data=zeros (trials, 4) ; 

12 

13 

%  Determine  initial  conditions 

14 

IC=  [200,  0,  4e— 7] ; 

15 

16 

%  Initialize  trials 

17 

for  trial=l : trials; 

18 

19 

20 

%  Generate  Coefficients 

21 

[ lambda, mu, k, A, p, c]  =  CoefCreator; 

22 

23 

%  Calculate  and  save  R  value 

24 

Data (trial, 4) = (k*lambda*p) / (c*A*mu) ; 

25 

26 

27 

%  Solve  the  ODE 

28 

%  The  output  is  the  matrix  with  time  and  solution.  The  imputs  are 

the 

29 

%  function  being  solved,  time  span,  initial  conditions,  options. 

and 

30 

%  other  inputs  to  the  function. 

31 

[t,y]  =  ode45 ( @Det_3CM, [ 0 , 100 ] , IC, [ ] , lambda, mu, k, A, p, c) ; 

32 

33 

34 

%  Determine  how  long  the  solution  vector  y  is 

35 

length_y=length (y) ; 

36 

37 

%  Take  the  end  values  of  the  solution  vector  to  y  and  store  them 

in  data 

38 

Data(trial,  [1:3] )  =y  ( length.y ,  : )  ; 

39 

40 

if  mod (trial, floor (trials /I 00) ) ==0 

41 

fprintf (' Percent  Completed:  %f%%,  Time  Elapsed:  %f  seconds  \n ' , 

trial/trials*100,  toe) 

42 

end 

43 

44 

end 
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Plot  ODE: 


1 

clear 

2 

clc 

3 

cl  f 

4 

5 

Data=ODE_Simulator ; 

6 

7 

hold  on 

8 

subplot (2 , 2 , 1 ) ,  scatter (Data ( : 

: , 4 ) , Data 

9 

Population ’ ) 

subplot (2, 2, 2) ,  scatter (Data ( : 

: ,  4 ) , Data 

10 

Population ' ) 

subplot (2, 2, 3) ,  scatter (Data ( : 

: ,  4 ) , Data 

Population ' ) 


1)  ) 

,  xlabel  ( 

'R 

Value ' ) , 

ylabel i 

( 'End 

T— Cell 

2)  ) 

,  xlabel  ( 

'R 

Value ' )  , 

ylabel i 

;  'End 

I— Cell 

3)  ) 

,  xlabel ( 

'R 

Value ' ) , 

ylabel i 

( 'End 

Virus 
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A. 0.3.  Discrete  Simulation  Code. 


1 

clear 

2 

clc 

3 

cl  f 

4 

5 

6 

7 

%Set  initial  values 

8 

T0=10000; 

9 

10=0; 

10 

V0=1 ; 

11 

I nit ia 1 Value s= [TO, 10, VO] ; 

12 

13 

14 

%Set  Coefficients 

15 

lambda  =  .272; 

16 

mu  =  .00136; 

17 

K  =  .00027; 

18 

A  =  .33; 

19 

p=100; 

20 

c=2 ; 

21 

22 

%Set  number  of  steps 

23 

steps=100000; 

24 

25 

%Set  total  time 

26 

Time=100 ; 

27 

28 

%Caclulate  time  interval 

29 

t=Time/steps; 

30 

31 

%Generate  "Storage"  Matrix  (will  be 

:  3x (Steps+1)  sized) 

32 

33 

Population=zeros  (steps+1, 3) ; 

34 

Population ( 1 , : ) =InitialValues; 

35 

36 

%Define  State  Changes 

37 

O 

O 

\ — i 

II 

i — i 
CO 

38 

S2=  [-1,  0,0]  ; 

39 

S3=  [-1, 1,0]  ; 

40 

o 

I — 1 

1 

o 

II 

CO 

41 

S5=  [ 0 ,  0,1]; 

42 

S6= [0, 0, —  1]  ; 

43 

CO 

o 

II 

o 

o 

o 

44 

45 

46 

%Define  number  of  trials; 

47 

trials=1000; 

48 

49 

%Set  storage  matricies  for  average 

value  and  end  states 

50 

Average=zeros (steps+1, 3) ; 

51 

EndValues=zeros (trials, 3) ; 

52 

53 

%Create  a  time  vector  for  plotting 

purposes 

54 

time= [ 0 : t : Time] ; 

55 

56 
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57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 


for  trial=l : trials ; 

Population ( 1 , : ) =InitialValues; 


for  k=l: steps 

%Calculate  probabilities  of  each  state  change 
pl=lambda*t ; 

p2=mu*Population (k, 1) *t; 
p3=K*Population (k, 1) ^Population (k, 3) *t ; 
p4=A*Population (k, 2) *t ; 
p5=p*Population (k, 2) *t; 
p6=c*Population (k, 3) *t; 
psum= (pl+p2+p3+p4+p5+p6) ; 
p0= (1— psum) *t; 


%Generate  a  uniform  random  number,  scale  it  to  the  sum  of  probabilities 


x=rand*psum; 


%State  change  test 
if  x<pl  &&  x  >  0 

Population (k+1 , : ) =Population (k, : ) +S1 ; 
elseif  x<pl+p2  &&  x  >  0 

Population (k+1, : ) =Population (k, : ) +S2 ; 
elseif  x<pl+p2+p3  &&  x  >  0 

Population (k+1, : ) =Population (k, : ) +S3; 
elseif  x<pl+p2+p3+p4  &&  x  >  0 

Population (k+1, : ) =Population (k, : ) +S4 ; 
elseif  X<pl+p2+p3+p4+p5  &&  x  >  0 

Population (k+1 , : ) =Population (k, : ) +S5 ; 
elseif  X<pl+p2+p3+p4+p5+p6  &&  x  >  0 

Population (k+1, : ) =Population (k, : ) +S6; 

else 

Population (k+1, : ) =Population (k, : ) +S0; 

end 


%Setting  the  floor 
%for  j=l:3 

%  if  Population (k+1, j ) <0 

%  Population (k+1, j ) =0; 

%  end 

%end 


end 


%Add  end  state  population  values  to  the  average  to  compute  the  average 
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118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135  end 

136 

137 

138  Average=Average/trials; 

139 

140  subplot  (2, 3, 1) ,  plot (time, Average (:, 1 ),' r ') ,  xlabel ( ' Time ' ) ,  ylabel (' Healthy  T— cells') 

141  subplot (2 , 3, 2 ) ,  plot (time, Average (:, 2 ),' r ') ,  xlabel (' Time ') ,  ylabel (' Infected  T— cells') 

142  subplot  (2, 3, 3) ,  plot (time, Average (:, 3) ,' r ') ,  xlabel (' Time ') ,  ylabel ( 'Virus  Particles') 

143  subplot  (2, 3, 4) ,  hist (EndValues (:, 1 ), 100,  title (' Healthy  T— Cell  Quantity')) 

144  subplot  (2, 3, 5) ,  hist (EndValues (:, 2 ), 100,  title (' Infected  T— Cell  Quantity')) 

145  subplot (2 ,  3, 6) ,  hist (EndValues (:, 3) , 100,  title ('Virus  Quantity')) 


%later 

Average=Average+Population; 


%Create  plots  of  T,I,  and  V  against  time 
hold  on 

subplot (2, 3,1),  plot (time, Population ( : , 1) ) 
hold  on 

subplot  (2, 3,2),  plot (time, Population ( : , 2) ) 
hold  on 

subplot  (2,3,3) ,  plot (time, Population ( : , 3) ) 
%Store  end  values  for  display  purposes 
EndValues (trial,  : ) =Population  (steps,  : ) ; 
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A. 0.4.  Stochastic  Code. 

Simulate  and  plot  ODE  with  random  valued  IV : 


1  function  StochasticSimulator 

2 

3  clear 

4  clc 

5  cl  f 

6 

7  %We  want  to  go  from  time  t=0  to  the  end  time,  t=T.  We  want  to  do  this  in  N 

8  %steps  of  length  L=T/N.  This  would  be  an  iterative  process.  We  start  with  the 

9  %populations  at  time  t=0 .  From  there,  we  can  calculate  populations  at  time 

10  %L  by  calculating  EV,  multiplying  it  by  L  and  taking  VAR,  multiplying  it  by 

11  %L,  then  multiplying  that  by  a  dW  of  unit  length.  We  would  start  with  a 

12  %matrix  with  three  columns  (for  T,  I,  &  V)  and  N+l  rows. 

13 

14 

15  %Coef f icients 

16  lambda  =  .272; 

17  mu  =  .00136; 

18  k  =  .00027; 

19  A  =  .33; 

20  p=100; 

21  c=2 ; 

22 

23  %Initial  Values  (T,I,V) 

24  InitialValues= [200, 0, 4e— 7] ; 

25 

26  R= (k*lambda*p) / (c*A*mu) 

27 

28  trials=2000; 

29  T=60;  %Total  Time  elapsed 

30  N=600;  %  number  of  trials 

31  dt=T/N; 

32  time=  [  0  :  dt :  T]  ; 

33  P=zeros (N+l,  3)  ; 

34  Average=zeros (N+l ,  3 )  ; 

35  EndValues=zeros (trials, 3) ; 

36  for  trial=l : trials; 

37  P  ( 1 ,  : ) =InitialValues ; 


38 


39 


40 


41 


for  i=l:N 


42 


44 


43 


k*P  (i,  1)  *p  (i,  3)  +A*P  (i,  2)  ,  0]  ;  [0, 0,p*P  (i, 2) +c*P  (i,  3)  ]  ]  ; 


45 


50 


46 


47 


48 


49 


51 


dW= [randn; randn; randn] ; 
dWt=dW*dt " ( . 5) ; 

%A=chol (VAR) ; 

A= (VAR) * . 5; 
dP=EV*dt+ (A) *dWt ; 

P ( i+1 ,  :  )=P  (i,  : )  +  [dP  (1)  ,dP  (2)  ,dP  (3)  ]  ; 


52 


53 


54 


%Absorbant  Boundrary  Conditions 
%P (i+1,  1)  =  max (  P (i  +  1,1),  eps)  ; 
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62 


55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 


%P (i+1,  2)  =  max (  P  (i  +  1,2),  0)  ; 
%P  (i+1,  3)  =  max (  P (i  +  1,3),  0)  ; 


%Ref lective 
P (i+1,  1)  = 
P (i+1,  2)  = 
P (i+1,  3)  = 


Boundrary  Conditions 
abs (  P (i  +  1, 1) )  ; 

abs (  P (i  +  1, 2) )  ; 

abs  (  P  (i  +  1, 3)  )  ; 


end 


%EV= [lambda— mu* P (l)-k*P (1) *P (3) ;k*P (1) *P (3)-A*P (2) ;p*P (2)-c*P (3) ] ; 
%VAR= [ [lambda— mu*P (l)-k*P (1) *P (3) ,  -k*P (1) *P (3) ,  0]  ;  [-k*P  (1) *P  (3) , 
k*P  (1)  *P  ( 3 )— A*P  (2) , 0] ;  [0, 0,p*P (2)-c*P (3) ] ] 


hold  on 

subplot  (2,3,1),  plot (time, P ( : , 1 ) ) 
hold  on 

subplot (2,3,2) ,  plot (time, P ( : , 2) ) 
hold  on 

subplot (2,3,3),  plot (time, P ( : , 3) ) 
Average=Average+P ; 

EndValues (trial, : ) =P (N+l, : ) ; 

end 


Average=Average/trials ; 

subplot (2, 3, 1) ,  plot (time, Average (:, 1 ),' r ') ,  xlabel ( ' Time ' ) ,  ylabel (' Healthy  T— cells') 
subplot (2, 3, 2) ,  plot (time, Average (:, 2 ),' r ') ,  xlabel (' Time ') ,  ylabel (' Infected  T— cells') 
subplot (2, 3, 3) ,  plot (time, Average (:, 3 ),' r ') ,  xlabel (' Time ') ,  ylabel (' Virus  Particles') 
subplot  (2, 3, 4) ,  hist (EndValues (:, 1 ), 100,  title (' Healthy  T— Cell  Quantity')) 
subplot (2, 3, 5) ,  hist (EndValues (:, 2 ), 100,  title (' Infected  T— Cell  Quantity')) 
subplot (2, 3, 6) ,  hist (EndValues (:, 3) , 100,  title ('Virus  Quantity')) 


